##Install Packages if Needed
if (!require("ggplot2")) install.packages("ggplot2")
if (!require("cowplot")) install.packages("cowplot")
if (!require("lubridate")) install.packages("lubridate")
if (!require("corrplot")) install.packages("corrplot")
if (!require("lme4")) install.packages("lme4")
if (!require("vegan")) install.packages("vegan")
if (!require("DHARMa")) install.packages("DHARMa")
if (!require("effectsize")) install.packages("effectsize")
if (!require("emmeans")) install.packages("emmeans")
if (!require("lmerTest")) install.packages("lmerTest")
if (!require("tidyr")) install.packages("tidyr")
if (!require("dplyr")) install.packages("dplyr")
if (!require("plotrix")) install.packages("plotrix")
if (!require("Rmisc")) install.packages("Rmisc")
if (!require("ggpubr")) install.packages("ggpubr")
##Load Packages
library(ggplot2) #Required for ggplots
library(cowplot) #Required for plotting panel figures
library(lubridate) #Required for date and time formatting
library(corrplot) #Required for correlation plot
library(lme4) #Required for mixed effects modeling
library(vegan) #Required for multivariate analysis PERM
library(DHARMa) #Required to check residuals of mixed effects models
library(effectsize) #Required for eta_squared effect sizes
library(emmeans) #Required for pairwise comparisons
library(lmerTest) #Required for p values with lme4 model summaries
library(tidyr) #Required for pivot_wider function
library(dplyr) #Required for group_by function
library(plotrix) #Required for Standard error function
library(Rmisc) #Required for SummarySE function
library(ggpubr) #Required for stat brackets
#Note: Run "Graphing Parameters" section from 01_ExperimentalSetup.Rmd file
##Growth
Length<-read.csv("Data/LengthFull.csv", header=TRUE)
Dates<-read.csv("Data/BonaireDates.csv", header=TRUE)
##Set date variables
Dates$Date<-as.Date(Dates$Date, format='%m/%d/%Y')
##Set factor variables
Length$Site<-factor(Length$Site, levels=c("KL", "SS"))
Length$Genotype<-factor(Length$Genotype, levels=c("AC8", "AC10", "AC12"))
Length$Orig<-factor(Length$Orig, levels=c("N", "T"))
Length$Origin<-factor(Length$Origin, levels=c("Native", "Transplant"))
##Add a Sample Set Variable
Length$Set<-paste(Length$Site, Length$Genotype, Length$Orig, sep=".")
Length$Set<-factor(Length$Set)
##Add Site.Orig variable
Length$Site.Orig<-paste(Length$Site, Length$Orig, sep=".")
Length$Site.Orig<-factor(Length$Site.Orig, levels=c("KL.N", "KL.T", "SS.N", "SS.T"))
##Bleaching Metrics
#Note: Physiological metrics calculated in 02_PhysiologyMetrics.R file
Thermal<-read.csv("Outputs/ThermData.csv", header=TRUE)
##Set factor variables
Thermal$TimeP<-factor(Thermal$TimeP, levels=c("IN", "TP1", "TP2", "TP3", "TP4"))
Thermal$Site<-factor(Thermal$Site, levels=c("KL", "SS"))
Thermal$Genotype<-factor(Thermal$Genotype, levels=c("AC8", "AC10", "AC12"))
Thermal$Orig<-factor(Thermal$Orig, levels=c("N", "T"))
Thermal$Origin<-factor(Thermal$Origin, levels=c("Native", "Transplant"))
##Subset Dates for TP 1 and TP 2 Growth
Growth.Dates_TP1.2<-subset(Dates, Task=="Growth" & TimeP=="TP1" |
Task=="Growth" & TimeP=="TP2" )
##Convert to Wide format
Growth.Dates_TP1.2<-Growth.Dates_TP1.2 %>% pivot_wider(names_from = TimeP, values_from = Date)
##Calculate Number of Days
#TP 1 to TP2
Growth.Dates_TP1.2$TP1.2_days<-time_length(Growth.Dates_TP1.2$TP2-Growth.Dates_TP1.2$TP1, unit="days")
##Merge Length Data with Growth Dates
#Merges by Site column
#Adds TP1.2_days column
Growth_TP1.2<-merge(Length, Growth.Dates_TP1.2[,c(-2)], all.x=TRUE)
##Calculate Linear Extension (cm)
Growth_TP1.2$Ext_cm<-Growth_TP1.2$TL.TP2_cm-Growth_TP1.2$TL.TP1_cm
##Calculate Growth_TP1.2 Rate (cm/day)
Growth_TP1.2$TLE_cm.day<-Growth_TP1.2$Ext_cm/Growth_TP1.2$TP1.2_days
##Remove Corals not Measured in T1 to T2
Growth_TP1.2<-Growth_TP1.2 %>% drop_na(TLE_cm.day)
##Sample Size per Set
Growth_TP1.2 %>%
dplyr::group_by(Set) %>%
dplyr::summarize(SampleSize = length(Set))
##Check for Outliers
ggplot(Growth_TP1.2, aes(x=Set, y=TLE_cm.day)) +
geom_boxplot(alpha=0.5, shape=2, outlier.shape = NA)+
geom_jitter(shape=16, position=position_jitter(0.1))+
theme(axis.text.x = element_text(angle = 90))
Potential outliers have been re-measured and therefore will be retained. n=10 per Site, Genotype, Origin group
#Check normality
hist(Growth_TP1.2$TLE_cm.day)
shapiro.test(Growth_TP1.2$TLE_cm.day)
Shapiro-Wilk normality test
data: Growth_TP1.2$TLE_cm.day
W = 0.92803, p-value = 7.262e-06
#Not normal
hist(log(Growth_TP1.2$TLE_cm.day+1))
shapiro.test(log(Growth_TP1.2$TLE_cm.day+1))
Shapiro-Wilk normality test
data: log(Growth_TP1.2$TLE_cm.day + 1)
W = 0.94335, p-value = 7.248e-05
##Still not normal
##Compare generalized linear mixed effects models
TLE_TP1.2.glmr.gaus<-glmer(TLE_cm.day~Origin*Site+(1|Genotype), data=Growth_TP1.2, family=gaussian(link="identity"))
Warning: calling glmer() with family=gaussian (identity link) as a shortcut to lmer() is deprecated; please call lmer() directly
TLE_TP1.2.glmr.gaus.l<-glmer(TLE_cm.day~Origin*Site+(1|Genotype), data=Growth_TP1.2, family=gaussian(link="log"))
TLE_TP1.2.glmr.gaus.i<-glmer(TLE_cm.day~Origin*Site+(1|Genotype), data=Growth_TP1.2, family=gaussian(link="inverse"))
TLE_TP1.2.glmr.gam.l<-glmer(TLE_cm.day~Origin*Site+(1|Genotype), data=Growth_TP1.2, family=Gamma(link="log"))
TLE_TP1.2.glmr.gam.i<-glmer(TLE_cm.day~Origin*Site+(1|Genotype), data=Growth_TP1.2, family=Gamma(link="inverse"))
TLE_TP1.2.glmr.inv<-glmer(TLE_cm.day~Origin*Site+(1|Genotype), data=Growth_TP1.2, family=inverse.gaussian(link="1/mu^2"))
boundary (singular) fit: see help('isSingular')
TLE_TP1.2.glmr.inv.l<-glmer(TLE_cm.day~Origin*Site+(1|Genotype), data=Growth_TP1.2, family=inverse.gaussian(link="log"))
TLE_TP1.2.glmr.inv.i<-glmer(TLE_cm.day~Origin*Site+(1|Genotype), data=Growth_TP1.2, family=inverse.gaussian(link="inverse"))
AIC(TLE_TP1.2.glmr.gaus, TLE_TP1.2.glmr.gaus.l, TLE_TP1.2.glmr.gaus.i, TLE_TP1.2.glmr.gam.l, TLE_TP1.2.glmr.gam.i, TLE_TP1.2.glmr.inv, TLE_TP1.2.glmr.inv.l, TLE_TP1.2.glmr.inv.i)
Gamma distribution with log-link has the lowest AIC.
##Check residuals
TLE_TP1.2.glmr.gam.l_res <- simulateResiduals(fittedModel = TLE_TP1.2.glmr.gam.l, plot = F)
plot(TLE_TP1.2.glmr.gam.l_res)
Compare residuals across models
##Check residuals
TLE_TP1.2.glmr.gaus.l_res <- simulateResiduals(fittedModel = TLE_TP1.2.glmr.gaus.l, plot = F)
plot(TLE_TP1.2.glmr.gaus.l_res)
Compare with log+1 transformed lmer model
##Model with log +1 transformation
#Function of Site and Origin, with Genotype as a Random effect
#Interaction between Site and Origin
TLE_TP1.2.lme<-lmer(log(TLE_cm.day+1)~Origin*Site+(1|Genotype), data=Growth_TP1.2)
##Check residuals
TLE_TP1.2.lme_res <- simulateResiduals(fittedModel = TLE_TP1.2.lme, plot = F)
plot(TLE_TP1.2.lme_res)
Better residuals with the log+1 transformed lmer model.
##Model results
summary(TLE_TP1.2.lme)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(TLE_cm.day + 1) ~ Origin * Site + (1 | Genotype)
Data: Growth_TP1.2
REML criterion at convergence: -233.2
Scaled residuals:
Min 1Q Median 3Q Max
-2.0026 -0.5991 -0.1118 0.5959 2.4282
Random effects:
Groups Name Variance Std.Dev.
Genotype (Intercept) 0.002965 0.05445
Residual 0.006632 0.08144
Number of obs: 120, groups: Genotype, 3
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.125179 0.034778 2.684719 3.599 0.04400 *
OriginTransplant -0.007538 0.021027 114.000000 -0.359 0.72063
SiteSS 0.062096 0.021027 114.000000 2.953 0.00382 **
OriginTransplant:SiteSS -0.035167 0.029737 114.000000 -1.183 0.23943
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) OrgnTr SiteSS
OrgnTrnspln -0.302
SiteSS -0.302 0.500
OrgnTrn:SSS 0.214 -0.707 -0.707
eta_squared(TLE_TP1.2.lme)
# Effect Size for ANOVA (Type III)
Parameter | Eta2 (partial) | 95% CI
-------------------------------------------
Origin | 0.02 | [0.00, 1.00]
Site | 0.07 | [0.01, 1.00]
Origin:Site | 0.01 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
TLE_TP1.2.res<-data.frame(summary(TLE_TP1.2.lme)$coefficients[-1,])
TLE_TP1.2.res$Predictor<-c("Origin", "Site", "Origin x Site")
TLE_TP1.2.res$EtaSq<-c(eta_squared(TLE_TP1.2.lme)$Eta2)
TLE_TP1.2.res$Response<-rep("Growth", nrow(TLE_TP1.2.res))
TLE_TP1.2.res$Timepoint<-rep("TP1_2", nrow(TLE_TP1.2.res))
##Pairwise Comparisons of Interest
#Native vs Transplant within a Site
TLE.pair_TP1.2<-emmeans(TLE_TP1.2.lme, pairwise~Origin | Site)
TLE.pair_TP1.2
$emmeans
Site = KL:
Origin emmean SE df lower.CL upper.CL
Native 0.125 0.0348 2.68 0.006770 0.244
Transplant 0.118 0.0348 2.68 -0.000768 0.236
Site = SS:
Origin emmean SE df lower.CL upper.CL
Native 0.187 0.0348 2.68 0.068866 0.306
Transplant 0.145 0.0348 2.68 0.026161 0.263
Degrees-of-freedom method: kenward-roger
Results are given on the log(mu + 1) (not the response) scale.
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
Native - Transplant 0.00754 0.021 114 0.359 0.7206
Site = SS:
contrast estimate SE df t.ratio p.value
Native - Transplant 0.04271 0.021 114 2.031 0.0446
Note: contrasts are still on the log(mu + 1) scale
Degrees-of-freedom method: kenward-roger
#Calculate standardized effect sizes for pairwise comparisons
TLE.pair.es_TP1.2<-data.frame(eff_size(TLE.pair_TP1.2, sigma=sigma(TLE_TP1.2.lme), edf=df.residual(TLE_TP1.2.lme)))
Since 'object' is a list, we are using the contrasts already present.
TLE.pair.es_TP1.2
TLE.pair.es_TP1.2<-TLE.pair.es_TP1.2 %>% dplyr::rename(contrast.se = contrast, SE.es = SE, df.es = df)
##Save Results
TLE.pair_TP1.2.res<-merge(data.frame(TLE.pair_TP1.2$contrasts), TLE.pair.es_TP1.2[,-c(1)])
TLE.pair_TP1.2.res$Response<-rep("Growth", nrow(TLE.pair_TP1.2.res))
TLE.pair_TP1.2.res$TimeP<-rep("TP1_2", nrow(TLE.pair_TP1.2.res))
##Add Significance levels
TLE.pair_TP1.2.res$Sig<-ifelse(TLE.pair_TP1.2.res$p.value<0.001, "***",
ifelse(TLE.pair_TP1.2.res$p.value<0.01, "**",
ifelse(TLE.pair_TP1.2.res$p.value<0.05, "*",
ifelse(TLE.pair_TP1.2.res$p.value<0.1, "-", NA))))
TLE.pair_TP1.2.res
##Summary statistics by Site and Origin
TP1.2_Growth.sum<-summarySE(Growth_TP1.2, measurevar="TLE_cm.day", groupvars=c("Site", "Origin", "Site.Orig"), na.rm=TRUE)
##Plot Average Growth Rate across Treatments
TP1.2_Growth.plot<-ggplot(TP1.2_Growth.sum, aes(x=Site, y=TLE_cm.day, colour=Site.Orig)) +
scale_colour_manual(values=Orig.colors.o)+
geom_errorbar(aes(ymin=TLE_cm.day-se, ymax=TLE_cm.day+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
geom_point(size=point.sz, position=position_dodge(width=0.5))+
theme_classic()+
theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"))+
labs(x="Site and Origin", y=expression(paste('Growth Rate (cm day'^-1*")")), colour="Origin")+
ylim(0, 0.25)+
annotate("text", x=2, y=0.24, label="*", size=sig.sz, fontface="bold")+
geom_bracket(xmin=1, xmax=2, y.position=0.01, label.size=levels.sz, label="**", inherit.aes = FALSE); TP1.2_Growth.plot
##Summary statistics by Genotype, Origin, and Site
TP1.2_Growth.sum<-summarySE(Growth_TP1.2, measurevar="TLE_cm.day", groupvars=c("Genotype", "Origin", "Site"), na.rm=TRUE)
##Calculate Percent Difference between Native and Transplant by Genotype
TP1.2_Growth.dif<-TP1.2_Growth.sum[which(TP1.2_Growth.sum$Origin=="Native"), c(1,3,5)]
TP1.2_Growth.dif<-TP1.2_Growth.dif %>% dplyr::rename(TLE_N = TLE_cm.day)
TP1.2_Growth.dif$TLE_T<-TP1.2_Growth.sum$TLE_cm.day[which(TP1.2_Growth.sum$Origin=="Transplant")]
##SS
TP1.2_Growth.dif_SS<-subset(TP1.2_Growth.dif, Site=="SS")
TP1.2_Growth.dif_SS$Perc.Inc<-((TP1.2_Growth.dif_SS$TLE_N-TP1.2_Growth.dif_SS$TLE_T)/TP1.2_Growth.dif_SS$TLE_T)*100
TP1.2_Growth.dif_SS
##Mean and Standard Error
mean(TP1.2_Growth.dif_SS$Perc.Inc)
[1] 47.19841
std.error(TP1.2_Growth.dif_SS$Perc.Inc)
[1] 25.89777
##KL
TP1.2_Growth.dif_KL<-subset(TP1.2_Growth.dif, Site=="KL")
TP1.2_Growth.dif_KL$Perc.Inc<-((TP1.2_Growth.dif_KL$TLE_N-TP1.2_Growth.dif_KL$TLE_T)/TP1.2_Growth.dif_KL$TLE_T)*100
TP1.2_Growth.dif_KL
##Mean and Standard Error
mean(TP1.2_Growth.dif_KL$Perc.Inc)
[1] 4.187086
std.error(TP1.2_Growth.dif_KL$Perc.Inc)
[1] 8.336196
##Subset Dates for TP 3 and TP 4 Growth
Growth.Dates_TP3.4<-subset(Dates, Task=="Growth" & TimeP=="TP3" |
Task=="Growth" & TimeP=="TP4" )
##Convert to Wide format
Growth.Dates_TP3.4<-Growth.Dates_TP3.4 %>% pivot_wider(names_from = TimeP, values_from = Date)
##Calculate Number of Days
#TP 3 to TP4
Growth.Dates_TP3.4$TP3.4_days<-time_length(Growth.Dates_TP3.4$TP4-Growth.Dates_TP3.4$TP3, unit="days")
##Merge Length Data with Growth Dates
#Merges by Site column
#Adds TP3.4_days column
Growth_TP3.4<-merge(Length, Growth.Dates_TP3.4[,c(-2)], all.x=TRUE)
##Calculate Linear Extension (cm)
Growth_TP3.4$Ext_cm<-Growth_TP3.4$TL.TP4_cm-Growth_TP3.4$TL.TP3_cm
##Calculate Growth_TP3.4 Rate (cm/day)
Growth_TP3.4$TLE_cm.day<-Growth_TP3.4$Ext_cm/Growth_TP3.4$TP3.4_days
##Remove Corals not Measured in T1 to T2
Growth_TP3.4<-Growth_TP3.4 %>% drop_na(TLE_cm.day)
##Sample Size per Set
Growth_TP3.4 %>%
dplyr::group_by(Set) %>%
dplyr::summarize(SampleSize = length(Set))
##Check for Outliers
ggplot(Growth_TP3.4, aes(x=Set, y=TLE_cm.day)) +
geom_boxplot(alpha=0.5, shape=2, outlier.shape = NA)+
geom_jitter(shape=16, position=position_jitter(0.1))+
theme(axis.text.x = element_text(angle = 90))
#Check normality
hist(Growth_TP3.4$TLE_cm.day)
shapiro.test(Growth_TP3.4$TLE_cm.day)
Shapiro-Wilk normality test
data: Growth_TP3.4$TLE_cm.day
W = 0.92767, p-value = 6.904e-06
#Not normal
hist(log(Growth_TP3.4$TLE_cm.day+1))
shapiro.test(log(Growth_TP3.4$TLE_cm.day+1))
Shapiro-Wilk normality test
data: log(Growth_TP3.4$TLE_cm.day + 1)
W = 0.94719, p-value = 0.0001351
##Still not normal
##Compare generalized linear mixed effects models
TLE_TP3.4.glmr.gaus<-glmer(TLE_cm.day~Origin*Site+(1|Genotype), data=Growth_TP3.4, family=gaussian(link="identity"))
Warning: calling glmer() with family=gaussian (identity link) as a shortcut to lmer() is deprecated; please call lmer() directly
TLE_TP3.4.glmr.gaus.l<-glmer(TLE_cm.day~Origin*Site+(1|Genotype), data=Growth_TP3.4, family=gaussian(link="log"))
TLE_TP3.4.glmr.gaus.i<-glmer(TLE_cm.day~Origin*Site+(1|Genotype), data=Growth_TP3.4, family=gaussian(link="inverse"))
TLE_TP3.4.glmr.gam.l<-glmer(TLE_cm.day~Origin*Site+(1|Genotype), data=Growth_TP3.4, family=Gamma(link="log"))
TLE_TP3.4.glmr.gam.i<-glmer(TLE_cm.day~Origin*Site+(1|Genotype), data=Growth_TP3.4, family=Gamma(link="inverse"))
TLE_TP3.4.glmr.inv<-glmer(TLE_cm.day~Origin*Site+(1|Genotype), data=Growth_TP3.4, family=inverse.gaussian(link="1/mu^2"))
TLE_TP3.4.glmr.inv.l<-glmer(TLE_cm.day~Origin*Site+(1|Genotype), data=Growth_TP3.4, family=inverse.gaussian(link="log"))
TLE_TP3.4.glmr.inv.i<-glmer(TLE_cm.day~Origin*Site+(1|Genotype), data=Growth_TP3.4, family=inverse.gaussian(link="inverse"))
AIC(TLE_TP3.4.glmr.gaus, TLE_TP3.4.glmr.gaus.l, TLE_TP3.4.glmr.gaus.i, TLE_TP3.4.glmr.gam.l, TLE_TP3.4.glmr.gam.i, TLE_TP3.4.glmr.inv, TLE_TP3.4.glmr.inv.l, TLE_TP3.4.glmr.inv.i)
Gamma distribution with log-link has the lowest AIC.
##Check residuals
TLE_TP3.4.glmr.gam.l_res <- simulateResiduals(fittedModel = TLE_TP3.4.glmr.gam.l, plot = F)
plot(TLE_TP3.4.glmr.gam.l_res)
Compare residuals across models
##Check residuals
TLE_TP3.4.glmr.gaus.l_res <- simulateResiduals(fittedModel = TLE_TP3.4.glmr.gaus.l, plot = F)
plot(TLE_TP3.4.glmr.gaus.l_res)
Compare with log+1 transformed lmer model
##Model with log +1 transformation
#Function of Site and Origin, with Genotype as a Random effect
#Interaction between Site and Origin
TLE_TP3.4.lme<-lmer(log(TLE_cm.day+1)~Origin*Site+(1|Genotype), data=Growth_TP3.4)
##Check residuals
TLE_TP3.4.lme_res <- simulateResiduals(fittedModel = TLE_TP3.4.lme, plot = F)
plot(TLE_TP3.4.lme_res)
Better residuals with the log+1 transformed lmer model.
##Model results
summary(TLE_TP3.4.lme)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(TLE_cm.day + 1) ~ Origin * Site + (1 | Genotype)
Data: Growth_TP3.4
REML criterion at convergence: -341.9
Scaled residuals:
Min 1Q Median 3Q Max
-2.27432 -0.60263 -0.01933 0.56722 2.45603
Random effects:
Groups Name Variance Std.Dev.
Genotype (Intercept) 0.002008 0.04481
Residual 0.002574 0.05074
Number of obs: 120, groups: Genotype, 3
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.16891 0.02748 2.38972 6.147 0.016262 *
OriginTransplant -0.02845 0.01310 114.00000 -2.172 0.031936 *
SiteSS -0.05868 0.01310 114.00000 -4.479 1.79e-05 ***
OriginTransplant:SiteSS 0.06284 0.01853 114.00000 3.392 0.000954 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) OrgnTr SiteSS
OrgnTrnspln -0.238
SiteSS -0.238 0.500
OrgnTrn:SSS 0.169 -0.707 -0.707
eta_squared(TLE_TP3.4.lme)
# Effect Size for ANOVA (Type III)
Parameter | Eta2 (partial) | 95% CI
-------------------------------------------
Origin | 9.01e-04 | [0.00, 1.00]
Site | 0.07 | [0.01, 1.00]
Origin:Site | 0.09 | [0.02, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
TLE_TP3.4.res<-data.frame(summary(TLE_TP3.4.lme)$coefficients[-1,])
TLE_TP3.4.res$Predictor<-c("Origin", "Site", "Origin x Site")
TLE_TP3.4.res$EtaSq<-c(eta_squared(TLE_TP3.4.lme)$Eta2)
TLE_TP3.4.res$Response<-rep("Growth", nrow(TLE_TP3.4.res))
TLE_TP3.4.res$Timepoint<-rep("TP3_4", nrow(TLE_TP3.4.res))
##Pairwise Comparisons of Interest
#Native vs Transplant within a Site
TLE.pair_TP3.4<-emmeans(TLE_TP3.4.lme, pairwise~Origin | Site)
TLE.pair_TP3.4
$emmeans
Site = KL:
Origin emmean SE df lower.CL upper.CL
Native 0.169 0.0275 2.39 0.06737 0.270
Transplant 0.140 0.0275 2.39 0.03892 0.242
Site = SS:
Origin emmean SE df lower.CL upper.CL
Native 0.110 0.0275 2.39 0.00869 0.212
Transplant 0.145 0.0275 2.39 0.04308 0.246
Degrees-of-freedom method: kenward-roger
Results are given on the log(mu + 1) (not the response) scale.
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
Native - Transplant 0.0285 0.0131 114 2.172 0.0319
Site = SS:
contrast estimate SE df t.ratio p.value
Native - Transplant -0.0344 0.0131 114 -2.625 0.0098
Note: contrasts are still on the log(mu + 1) scale
Degrees-of-freedom method: kenward-roger
#Calculate standardized effect sizes for pairwise comparisons
TLE.pair.es_TP3.4<-data.frame(eff_size(TLE.pair_TP3.4, sigma=sigma(TLE_TP3.4.lme), edf=df.residual(TLE_TP3.4.lme)))
Since 'object' is a list, we are using the contrasts already present.
TLE.pair.es_TP3.4
TLE.pair.es_TP3.4<-TLE.pair.es_TP3.4 %>% dplyr::rename(contrast.se = contrast, SE.es = SE, df.es = df)
##Save Results
TLE.pair_TP3.4.res<-merge(data.frame(TLE.pair_TP3.4$contrasts), TLE.pair.es_TP3.4[,-c(1)])
TLE.pair_TP3.4.res$Response<-rep("Growth", nrow(TLE.pair_TP3.4.res))
TLE.pair_TP3.4.res$TimeP<-rep("TP3_4", nrow(TLE.pair_TP3.4.res))
##Add Significance levels
TLE.pair_TP3.4.res$Sig<-ifelse(TLE.pair_TP3.4.res$p.value<0.001, "***",
ifelse(TLE.pair_TP3.4.res$p.value<0.01, "**",
ifelse(TLE.pair_TP3.4.res$p.value<0.05, "*",
ifelse(TLE.pair_TP3.4.res$p.value<0.1, "-", NA))))
TLE.pair_TP3.4.res
##Summary statistics by Site and Origin
TP3.4_Growth.sum<-summarySE(Growth_TP3.4, measurevar="TLE_cm.day", groupvars=c("Site", "Origin", "Site.Orig"), na.rm=TRUE)
##Plot Average Growth Rate across Treatments
TP3.4_Growth.plot<-ggplot(TP3.4_Growth.sum, aes(x=Site, y=TLE_cm.day, colour=Site.Orig)) +
scale_colour_manual(values=Orig.colors.o)+
geom_errorbar(aes(ymin=TLE_cm.day-se, ymax=TLE_cm.day+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
geom_point(size=point.sz, position=position_dodge(width=0.5))+
theme_classic()+
theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"))+
labs(x="Site and Origin", y=expression(paste('Growth Rate (cm day'^-1*")")), colour="Origin")+
ylim(0, 0.25)+
annotate("text", x=1, y=0.22, label="*", size=sig.sz, fontface="bold")+
annotate("text", x=2, y=0.18, label="**", size=sig.sz, fontface="bold")+
geom_bracket(xmin=1, xmax=2, y.position=0.01, label.size=levels.sz, label="***", inherit.aes = FALSE); TP3.4_Growth.plot
##Summary statistics by Genotype, Origin, and Site
TP3.4_Growth.sum<-summarySE(Growth_TP3.4, measurevar="TLE_cm.day", groupvars=c("Genotype", "Origin", "Site"), na.rm=TRUE)
##Calculate Percent Difference between Native and Transplant by Genotype
TP3.4_Growth.dif<-TP3.4_Growth.sum[which(TP3.4_Growth.sum$Origin=="Native"), c(1,3,5)]
TP3.4_Growth.dif<-TP3.4_Growth.dif %>% dplyr::rename(TLE_N = TLE_cm.day)
TP3.4_Growth.dif$TLE_T<-TP3.4_Growth.sum$TLE_cm.day[which(TP3.4_Growth.sum$Origin=="Transplant")]
##SS
TP3.4_Growth.dif_SS<-subset(TP3.4_Growth.dif, Site=="SS")
TP3.4_Growth.dif_SS$Perc.Inc<-((TP3.4_Growth.dif_SS$TLE_T-TP3.4_Growth.dif_SS$TLE_N)/TP3.4_Growth.dif_SS$TLE_N)*100
TP3.4_Growth.dif_SS
##Mean and Standard Error
mean(TP3.4_Growth.dif_SS$Perc.Inc)
[1] 34.22395
std.error(TP3.4_Growth.dif_SS$Perc.Inc)
[1] 18.61694
##KL
TP3.4_Growth.dif_KL<-subset(TP3.4_Growth.dif, Site=="KL")
TP3.4_Growth.dif_KL$Perc.Inc<-((TP3.4_Growth.dif_KL$TLE_N-TP3.4_Growth.dif_KL$TLE_T)/TP3.4_Growth.dif_KL$TLE_T)*100
TP3.4_Growth.dif_KL
##Mean and Standard Error
mean(TP3.4_Growth.dif_KL$Perc.Inc)
[1] 21.37099
std.error(TP3.4_Growth.dif_KL$Perc.Inc)
[1] 4.101183
##Growth Data dataframe
GrowthData<-Growth_TP1.2
##Specify Timepoints
GrowthData<-GrowthData %>% dplyr::rename(TP1.2_Ext_cm = Ext_cm, TP1.2_TLE_cm.day = TLE_cm.day)
##Merge with Timepoints 3 to 4
GrowthData<-merge(GrowthData, Growth_TP3.4, all=TRUE)
##Specify Timepoints
GrowthData<-GrowthData %>% dplyr::rename(TP3.4_Ext_cm = Ext_cm, TP3.4_TLE_cm.day = TLE_cm.day)
##Growth Data
write.csv(GrowthData, "Outputs/GrowthData.csv", row.names=FALSE)
##Control
Therm_C<-subset(Thermal, Treat=="C")
##Heated
Therm_H<-subset(Thermal, Treat=="H")
Calculating retention as proportion remaining relative to corresponding control levels (0-1) for each site and genotype at each timepoint.
##Calculate averages for Control Treatment for each Site, Genotype, and Timepoint
Therm_C<-na.omit(Therm_C)
names(Therm_C)
[1] "ID" "RandN" "TimeP" "Site" "Genotype" "Treat"
[7] "Treatment" "Orig" "Origin" "Set" "Site.Orig" "SA_cm2"
[13] "Chl_ug.cm2" "Sym10.6_cm2" "Fv_Fm"
Therm_C.a<-aggregate(Therm_C[,c(13:15)], list(Therm_C$Site, Therm_C$Genotype, Therm_C$TimeP, Therm_C$Origin), mean, na.action = na.omit)
names(Therm_C.a)[1:4]<-c("Site", "Genotype", "TimeP", "Origin")
names(Therm_C.a)[5:7]<-paste(names(Therm_C.a)[5:7], "C", sep="_")
##Merge Control Averages with Heated Samples
#Merges by Site, Genotype, and Timepoint
Therm_H<-merge(Therm_H, Therm_C.a, all.x=TRUE )
##Calculate Proportion Retained for each Bleaching Metric relative to Control
Therm_H$Chl.prop<-round((Therm_H$Chl_ug.cm2/Therm_H$Chl_ug.cm2_C), 4)
Therm_H$Sym.prop<-round((Therm_H$Sym10.6_cm2/Therm_H$Sym10.6_cm2_C), 4)
Therm_H$PAM.prop<-round((Therm_H$Fv_Fm/Therm_H$Fv_Fm_C), 4)
##Set values >1 to 1
Therm_H$Chl.prop[which(Therm_H$Chl.prop>1)]<-1.0000
Therm_H$Sym.prop[which(Therm_H$Sym.prop>1)]<-1.0000
Therm_H$PAM.prop[which(Therm_H$PAM.prop>1)]<-1.0000
##Subset Initial Timepoint
Thermal_IN<-subset(Thermal, TimeP=="IN")
Thermal_IN$Site.Treat<-paste(Thermal_IN$Site, Thermal_IN$Treat, sep=".")
Thermal_IN$Site.Treat<-factor(Thermal_IN$Site.Treat, levels=c("KL.C", "KL.H", "SS.C", "SS.H"))
Therm_IN<-subset(Therm_H, TimeP=="IN")
#Check normality
hist(Thermal_IN$Sym10.6_cm2)
shapiro.test(Thermal_IN$Sym10.6_cm2)
Shapiro-Wilk normality test
data: Thermal_IN$Sym10.6_cm2
W = 0.90824, p-value = 0.001334
#Not normal
hist(log(Thermal_IN$Sym10.6_cm2+1))
shapiro.test(log(Thermal_IN$Sym10.6_cm2+1))
Shapiro-Wilk normality test
data: log(Thermal_IN$Sym10.6_cm2 + 1)
W = 0.95667, p-value = 0.07936
#Normal
##Model
#Function of Site and Treatment, with Genotype as a Random effect
Sym.lme_IN<-lmer(log(Sym10.6_cm2+1)~Site*Treatment+(1|Genotype), data=Thermal_IN)
##Check residuals
Sym.lme_res_IN <- simulateResiduals(fittedModel = Sym.lme_IN, plot = F)
plot(Sym.lme_res_IN)
##Model results
summary(Sym.lme_IN)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(Sym10.6_cm2 + 1) ~ Site * Treatment + (1 | Genotype)
Data: Thermal_IN
REML criterion at convergence: -29.3
Scaled residuals:
Min 1Q Median 3Q Max
-1.79325 -0.66013 0.08255 0.57806 2.12486
Random effects:
Groups Name Variance Std.Dev.
Genotype (Intercept) 0.0008222 0.02867
Residual 0.0230525 0.15183
Number of obs: 47, groups: Genotype, 3
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.38823 0.04685 14.06727 8.286 8.77e-07 ***
SiteSS 0.19508 0.06198 41.00308 3.147 0.00307 **
TreatmentHeat -0.20370 0.06340 41.10985 -3.213 0.00256 **
SiteSS:TreatmentHeat -0.01724 0.08867 41.05799 -0.194 0.84677
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) SiteSS TrtmnH
SiteSS -0.661
TreatmentHt -0.647 0.489
StSS:TrtmnH 0.462 -0.699 -0.715
eta_squared(Sym.lme_IN)
# Effect Size for ANOVA (Type III)
Parameter | Eta2 (partial) | 95% CI
----------------------------------------------
Site | 0.30 | [0.12, 1.00]
Treatment | 0.36 | [0.17, 1.00]
Site:Treatment | 9.20e-04 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
emmeans(Sym.lme_IN, pairwise ~ Site | Treatment)
$emmeans
Treatment = Control:
Site emmean SE df lower.CL upper.CL
KL 0.388 0.0469 14.1 0.2878 0.489
SS 0.583 0.0469 14.1 0.4829 0.684
Treatment = Heat:
Site emmean SE df lower.CL upper.CL
KL 0.185 0.0488 15.6 0.0808 0.288
SS 0.362 0.0469 14.1 0.2619 0.463
Degrees-of-freedom method: kenward-roger
Results are given on the log(mu + 1) (not the response) scale.
Confidence level used: 0.95
$contrasts
Treatment = Control:
contrast estimate SE df t.ratio p.value
KL - SS -0.195 0.0620 41.0 -3.147 0.0031
Treatment = Heat:
contrast estimate SE df t.ratio p.value
KL - SS -0.178 0.0635 41.1 -2.801 0.0077
Note: contrasts are still on the log(mu + 1) scale
Degrees-of-freedom method: kenward-roger
emmeans(Sym.lme_IN, pairwise ~ Treatment | Site)
$emmeans
Site = KL:
Treatment emmean SE df lower.CL upper.CL
Control 0.388 0.0469 14.1 0.2878 0.489
Heat 0.185 0.0488 15.6 0.0808 0.288
Site = SS:
Treatment emmean SE df lower.CL upper.CL
Control 0.583 0.0469 14.1 0.4829 0.684
Heat 0.362 0.0469 14.1 0.2619 0.463
Degrees-of-freedom method: kenward-roger
Results are given on the log(mu + 1) (not the response) scale.
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
Control - Heat 0.204 0.0635 41.1 3.209 0.0026
Site = SS:
contrast estimate SE df t.ratio p.value
Control - Heat 0.221 0.0620 41.0 3.564 0.0009
Note: contrasts are still on the log(mu + 1) scale
Degrees-of-freedom method: kenward-roger
##Save model results
Sym_IN.res<-data.frame(summary(Sym.lme_IN)$coefficients[-1,])
##Add pairwise
names(Sym_IN.res)<-names(data.frame(emmeans(Sym.lme_IN, pairwise ~ Site | Treatment)$contrasts))[-c(1:2)]
Sym_IN.res<-rbind(Sym_IN.res,
data.frame(rbind(emmeans(Sym.lme_IN, pairwise ~ Site | Treatment)$contrasts[1], emmeans(Sym.lme_IN, pairwise ~ Site | Treatment)$contrasts[2], emmeans(Sym.lme_IN, pairwise ~ Treatment | Site)$contrasts[1], emmeans(Sym.lme_IN, pairwise ~ Treatment | Site)$contrasts[2]))[,-c(1:3)])
##Metadata
Sym_IN.res$Predictor<-c("Site", "Treatment", "Site x Treatment", "Control Site", "Heated Site", "KL Treatment", "SS Treatment")
Sym_IN.res$Response<-rep("Symbionts", nrow(Sym_IN.res))
Sym_IN.res$EtaSq<-c(eta_squared(Sym.lme_IN)$Eta2, rep(NA, 4))
##Summary statistics by Site
IN_Sym.sum<-summarySE(Thermal_IN, measurevar="Sym10.6_cm2", groupvars=c("Site", "Treatment", "Site.Treat"), na.rm=TRUE)
##Plot Average Symbionts across Treatments
IN_Sym.plot<-ggplot(IN_Sym.sum, aes(x=Site.Treat, y=Sym10.6_cm2, colour=Site)) +
scale_colour_manual(values=Site.colors.o)+
geom_errorbar(aes(ymin=Sym10.6_cm2-se, ymax=Sym10.6_cm2+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
geom_point(aes(shape=Treatment), size=point.sz, position=position_dodge(width=0.5))+
scale_shape_manual(values=c(1,16))+
theme_classic()+
theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"))+
labs(x="Site and Treatment", y=expression(paste('Symbiont Density ('*10^6,'cells cm'^-2*")")), colour="Site")+
ylim(0, 1)+
annotate("text", x=.6, y=1, hjust=0, label="Site **", size=sig.sz, fontface="bold")+
annotate("text", x=.6, y=.94, hjust=0, label="Treatment **", size=sig.sz, fontface="bold"); IN_Sym.plot
#Check normality
hist(Thermal_IN$Chl_ug.cm2)
shapiro.test(Thermal_IN$Chl_ug.cm2)
Shapiro-Wilk normality test
data: Thermal_IN$Chl_ug.cm2
W = 0.88059, p-value = 0.0002132
#Not normal
hist(log(Thermal_IN$Chl_ug.cm2+1))
shapiro.test(log(Thermal_IN$Chl_ug.cm2+1))
Shapiro-Wilk normality test
data: log(Thermal_IN$Chl_ug.cm2 + 1)
W = 0.92511, p-value = 0.005653
#Still not normal but less skewed
##Model
#Function of Site and Treatment, with Genotype as a Random effect
Chl.lme_IN<-lmer(log(Chl_ug.cm2+1)~Site*Treatment+(1|Genotype), data=Thermal_IN)
##Check residuals
Chl.lme_res_IN <- simulateResiduals(fittedModel = Chl.lme_IN, plot = F)
plot(Chl.lme_res_IN)
##Model results
summary(Chl.lme_IN)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(Chl_ug.cm2 + 1) ~ Site * Treatment + (1 | Genotype)
Data: Thermal_IN
REML criterion at convergence: -69.1
Scaled residuals:
Min 1Q Median 3Q Max
-2.50589 -0.51517 0.07352 0.50169 2.35158
Random effects:
Groups Name Variance Std.Dev.
Genotype (Intercept) 0.006137 0.07834
Residual 0.007930 0.08905
Number of obs: 46, groups: Genotype, 3
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.45648 0.05202 2.95779 8.774 0.00329 **
SiteSS 0.27810 0.03635 39.99255 7.650 2.38e-09 ***
TreatmentHeat -0.37295 0.03721 40.00578 -10.024 1.80e-12 ***
SiteSS:TreatmentHeat -0.16117 0.05265 40.01278 -3.061 0.00393 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) SiteSS TrtmnH
SiteSS -0.349
TreatmentHt -0.341 0.489
StSS:TrtmnH 0.241 -0.691 -0.707
eta_squared(Chl.lme_IN)
# Effect Size for ANOVA (Type III)
Parameter | Eta2 (partial) | 95% CI
----------------------------------------------
Site | 0.58 | [0.41, 1.00]
Treatment | 0.88 | [0.82, 1.00]
Site:Treatment | 0.19 | [0.04, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
emmeans(Chl.lme_IN, pairwise ~ Site | Treatment)
$emmeans
Treatment = Control:
Site emmean SE df lower.CL upper.CL
KL 0.4565 0.0520 2.97 0.2899 0.623
SS 0.7346 0.0520 2.97 0.5680 0.901
Treatment = Heat:
Site emmean SE df lower.CL upper.CL
KL 0.0835 0.0526 3.11 -0.0808 0.248
SS 0.2005 0.0526 3.11 0.0361 0.365
Degrees-of-freedom method: kenward-roger
Results are given on the log(mu + 1) (not the response) scale.
Confidence level used: 0.95
$contrasts
Treatment = Control:
contrast estimate SE df t.ratio p.value
KL - SS -0.278 0.0364 40 -7.650 <.0001
Treatment = Heat:
contrast estimate SE df t.ratio p.value
KL - SS -0.117 0.0381 40 -3.070 0.0038
Note: contrasts are still on the log(mu + 1) scale
Degrees-of-freedom method: kenward-roger
emmeans(Chl.lme_IN, pairwise ~ Treatment | Site)
$emmeans
Site = KL:
Treatment emmean SE df lower.CL upper.CL
Control 0.4565 0.0520 2.97 0.2899 0.623
Heat 0.0835 0.0526 3.11 -0.0808 0.248
Site = SS:
Treatment emmean SE df lower.CL upper.CL
Control 0.7346 0.0520 2.97 0.5680 0.901
Heat 0.2005 0.0526 3.11 0.0361 0.365
Degrees-of-freedom method: kenward-roger
Results are given on the log(mu + 1) (not the response) scale.
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
Control - Heat 0.373 0.0372 40 10.022 <.0001
Site = SS:
contrast estimate SE df t.ratio p.value
Control - Heat 0.534 0.0372 40 14.353 <.0001
Note: contrasts are still on the log(mu + 1) scale
Degrees-of-freedom method: kenward-roger
##Save model results
Chl_IN.res<-data.frame(summary(Chl.lme_IN)$coefficients[-1,])
##Add pairwise
names(Chl_IN.res)<-names(data.frame(emmeans(Chl.lme_IN, pairwise ~ Site | Treatment)$contrasts))[-c(1:2)]
Chl_IN.res<-rbind(Chl_IN.res,
data.frame(rbind(emmeans(Chl.lme_IN, pairwise ~ Site | Treatment)$contrasts[1], emmeans(Chl.lme_IN, pairwise ~ Site | Treatment)$contrasts[2], emmeans(Chl.lme_IN, pairwise ~ Treatment | Site)$contrasts[1], emmeans(Chl.lme_IN, pairwise ~ Treatment | Site)$contrasts[2]))[,-c(1:3)])
##Metadata
Chl_IN.res$Predictor<-c("Site", "Treatment", "Site x Treatment", "Control Site", "Heated Site", "KL Treatment", "SS Treatment")
Chl_IN.res$Response<-rep("Chlorophyll", nrow(Chl_IN.res))
Chl_IN.res$EtaSq<-c(eta_squared(Chl.lme_IN)$Eta2, rep(NA, 4))
##Summary statistics by Site
IN_Chl.sum<-summarySE(Thermal_IN, measurevar="Chl_ug.cm2", groupvars=c("Site", "Treatment", "Site.Treat"), na.rm=TRUE)
##Plot Average Chlorophyll across Treatments
IN_Chl.plot<-ggplot(IN_Chl.sum, aes(x=Site.Treat, y=Chl_ug.cm2, colour=Site)) +
scale_colour_manual(values=Site.colors.o)+
geom_errorbar(aes(ymin=Chl_ug.cm2-se, ymax=Chl_ug.cm2+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
geom_point(aes(shape=Treatment), size=point.sz, position=position_dodge(width=0.5))+
scale_shape_manual(values=c(1,16))+
theme_classic()+
theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"))+
labs(x="Site and Treatment", y=expression(paste('Total Chlorophyll (\u03BCg cm'^-2*")")), colour="Site")+
ylim(0, 1.25)+
annotate("text", x=.6, y=1.2, hjust=0, label="Site ***", size=sig.sz, fontface="bold")+
annotate("text", x=.6, y=1.12, hjust=0, label="Treatment ***", size=sig.sz, fontface="bold")+
annotate("text", x=.6, y=1.04, hjust=0, label="Site x Treatment **", size=sig.sz, fontface="bold"); IN_Chl.plot
#Check normality
hist(Thermal_IN$Fv_Fm)
shapiro.test(Thermal_IN$Fv_Fm)
Shapiro-Wilk normality test
data: Thermal_IN$Fv_Fm
W = 0.87006, p-value = 9.097e-05
#Not normal
hist(Thermal_IN$Fv_Fm^2)
shapiro.test(Thermal_IN$Fv_Fm^2)
Shapiro-Wilk normality test
data: Thermal_IN$Fv_Fm^2
W = 0.90968, p-value = 0.001489
#Still not normal but less skewed
##Model
#Function of Site and Treatment, with Genotype as a Random effect
PAM.lme_IN<-lmer(Fv_Fm^2~Site*Treatment+(1|Genotype), data=Thermal_IN)
boundary (singular) fit: see help('isSingular')
##Check residuals
PAM.lme_res_IN <- simulateResiduals(fittedModel = PAM.lme_IN, plot = F)
plot(PAM.lme_res_IN)
##Model results
summary(PAM.lme_IN)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Fv_Fm^2 ~ Site * Treatment + (1 | Genotype)
Data: Thermal_IN
REML criterion at convergence: -156.1
Scaled residuals:
Min 1Q Median 3Q Max
-3.2454 -0.5025 0.1248 0.7046 1.5447
Random effects:
Groups Name Variance Std.Dev.
Genotype (Intercept) 8.057e-21 8.976e-11
Residual 1.235e-03 3.514e-02
Number of obs: 47, groups: Genotype, 3
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.38027 0.01015 43.00000 37.482 < 2e-16 ***
SiteSS 0.02899 0.01435 43.00000 2.021 0.0496 *
TreatmentHeat -0.07814 0.01467 43.00000 -5.327 3.45e-06 ***
SiteSS:TreatmentHeat 0.03376 0.02052 43.00000 1.645 0.1072
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) SiteSS TrtmnH
SiteSS -0.707
TreatmentHt -0.692 0.489
StSS:TrtmnH 0.494 -0.699 -0.715
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
eta_squared(PAM.lme_IN)
# Effect Size for ANOVA (Type III)
Parameter | Eta2 (partial) | 95% CI
----------------------------------------------
Site | 0.32 | [0.14, 1.00]
Treatment | 0.45 | [0.27, 1.00]
Site:Treatment | 0.06 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
emmeans(PAM.lme_IN, pairwise ~ Site | Treatment)
$emmeans
Treatment = Control:
Site emmean SE df lower.CL upper.CL
KL 0.380 0.0101 21.6 0.359 0.401
SS 0.409 0.0101 21.6 0.388 0.430
Treatment = Heat:
Site emmean SE df lower.CL upper.CL
KL 0.302 0.0106 23.5 0.280 0.324
SS 0.365 0.0101 21.6 0.344 0.386
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
$contrasts
Treatment = Control:
contrast estimate SE df t.ratio p.value
KL - SS -0.0290 0.0143 41.0 -2.021 0.0499
Treatment = Heat:
contrast estimate SE df t.ratio p.value
KL - SS -0.0627 0.0147 41.2 -4.268 0.0001
Degrees-of-freedom method: kenward-roger
emmeans(PAM.lme_IN, pairwise ~ Treatment | Site)
$emmeans
Site = KL:
Treatment emmean SE df lower.CL upper.CL
Control 0.380 0.0101 21.6 0.359 0.401
Heat 0.302 0.0106 23.5 0.280 0.324
Site = SS:
Treatment emmean SE df lower.CL upper.CL
Control 0.409 0.0101 21.6 0.388 0.430
Heat 0.365 0.0101 21.6 0.344 0.386
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
Control - Heat 0.0781 0.0147 41.2 5.316 <.0001
Site = SS:
contrast estimate SE df t.ratio p.value
Control - Heat 0.0444 0.0143 41.0 3.094 0.0036
Degrees-of-freedom method: kenward-roger
##Save model results
PAM_IN.res<-data.frame(summary(PAM.lme_IN)$coefficients[-1,])
##Add pairwise
names(PAM_IN.res)<-names(data.frame(emmeans(PAM.lme_IN, pairwise ~ Site | Treatment)$contrasts))[-c(1:2)]
PAM_IN.res<-rbind(PAM_IN.res,
data.frame(rbind(emmeans(PAM.lme_IN, pairwise ~ Site | Treatment)$contrasts[1], emmeans(PAM.lme_IN, pairwise ~ Site | Treatment)$contrasts[2], emmeans(PAM.lme_IN, pairwise ~ Treatment | Site)$contrasts[1], emmeans(PAM.lme_IN, pairwise ~ Treatment | Site)$contrasts[2]))[,-c(1:3)])
##Metadata
PAM_IN.res$Predictor<-c("Site", "Treatment", "Site x Treatment", "Control Site", "Heated Site", "KL Treatment", "SS Treatment")
PAM_IN.res$Response<-rep("FvFm", nrow(PAM_IN.res))
PAM_IN.res$EtaSq<-c(eta_squared(PAM.lme_IN)$Eta2, rep(NA, 4))
Not normal, but model residuals are OK.
##Summary statistics by Site
IN_PAM.sum<-summarySE(Thermal_IN, measurevar="Fv_Fm", groupvars=c("Site", "Treatment", "Site.Treat"), na.rm=TRUE)
##Plot Average FvFm across Treatments
IN_PAM.plot<-ggplot(IN_PAM.sum, aes(x=Site.Treat, y=Fv_Fm, colour=Site)) +
scale_colour_manual(values=Site.colors.o)+
geom_errorbar(aes(ymin=Fv_Fm-se, ymax=Fv_Fm+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
geom_point(aes(shape=Treatment), size=point.sz, position=position_dodge(width=0.5))+
scale_shape_manual(values=c(1,16))+
theme_classic()+
theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"))+
labs(x="Site and Treatment", y="Photochemical Efficiency (Fv/Fm)", colour="Site")+
ylim(0.5, 0.655)+
annotate("text", x=.6, y=.655, hjust=0, label="Site *", size=sig.sz, fontface="bold")+
annotate("text", x=.6, y=.645, hjust=0, label="Treatment ***", size=sig.sz, fontface="bold"); IN_PAM.plot
#Check normality
hist(Therm_IN$Sym.prop)
shapiro.test(Therm_IN$Sym.prop)
Shapiro-Wilk normality test
data: Therm_IN$Sym.prop
W = 0.93042, p-value = 0.1117
#Normal
##Model
#Function of Site, with Genotype as a Random effect
Tol_Sym.lme_IN<-lmer(Sym.prop~Site+(1|Genotype), data=Therm_IN)
##Check residuals
Tol_Sym.lme_res_IN <- simulateResiduals(fittedModel = Tol_Sym.lme_IN, plot = F)
plot(Tol_Sym.lme_res_IN)
##Model results
summary(Tol_Sym.lme_IN)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Sym.prop ~ Site + (1 | Genotype)
Data: Therm_IN
REML criterion at convergence: -2.8
Scaled residuals:
Min 1Q Median 3Q Max
-1.6802 -0.6578 -0.1521 0.6805 1.9852
Random effects:
Groups Name Variance Std.Dev.
Genotype (Intercept) 0.07818 0.2796
Residual 0.03039 0.1743
Number of obs: 23, groups: Genotype, 3
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.46139 0.16983 2.21175 2.717 0.1014
SiteSS 0.14316 0.07291 19.00328 1.964 0.0644 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr)
SiteSS -0.225
eta_squared(Tol_Sym.lme_IN)
# Effect Size for ANOVA (Type III)
Parameter | Eta2 (partial) | 95% CI
-----------------------------------------
Site | 0.17 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_Sym_IN.res<-data.frame(summary(Tol_Sym.lme_IN)$coefficients)
names(Tol_Sym_IN.res)<-names(Sym_IN.res)[1:5]
Tol_Sym_IN.res$Predictor<-rep("Site", nrow(Tol_Sym_IN.res))
Tol_Sym_IN.res$Response<-rep("Symbiont Retention", nrow(Tol_Sym_IN.res))
Tol_Sym_IN.res$EtaSq<-c(eta_squared(Tol_Sym.lme_IN)$Eta2)
##Combine results
Sym_IN.res<-rbind(Sym_IN.res, Tol_Sym_IN.res[2,])
##Summary statistics by Site
IN_TolSym.sum<-summarySE(Therm_IN, measurevar="Sym.prop", groupvars=c("Site"), na.rm=TRUE)
##Plot Average Symbiont Retention across Treatments
IN_TolSym.plot<-ggplot(IN_TolSym.sum, aes(x=Site, y=Sym.prop, colour=Site)) +
scale_colour_manual(values=Site.colors.o)+
geom_errorbar(aes(ymin=Sym.prop-se, ymax=Sym.prop+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
geom_point(size=point.sz, position=position_dodge(width=0.5))+
theme_classic()+
theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"))+
labs(x="Site", y="Symbiont Retention", colour="Site")+
ylim(0, 1)+
annotate("text", x=1.5, y=0.75, label="-", size=sig.sz, fontface="bold"); IN_TolSym.plot
#Check normality
hist(Therm_IN$Chl.prop)
shapiro.test(Therm_IN$Chl.prop)
Shapiro-Wilk normality test
data: Therm_IN$Chl.prop
W = 0.98432, p-value = 0.9694
#Normal
##Model
#Function of Site, with Genotype as a Random effect
Tol_Chl.lme_IN<-lmer(Chl.prop~Site+(1|Genotype), data=Therm_IN)
##Check residuals
Tol_Chl.lme_res_IN <- simulateResiduals(fittedModel = Tol_Chl.lme_IN, plot = F)
plot(Tol_Chl.lme_res_IN)
##Model results
summary(Tol_Chl.lme_IN)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Chl.prop ~ Site + (1 | Genotype)
Data: Therm_IN
REML criterion at convergence: -47.2
Scaled residuals:
Min 1Q Median 3Q Max
-1.3909 -0.5550 -0.1823 0.5284 1.9947
Random effects:
Groups Name Variance Std.Dev.
Genotype (Intercept) 0.004576 0.06765
Residual 0.003436 0.05861
Number of obs: 22, groups: Genotype, 3
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.16552 0.04290 2.40867 3.858 0.045 *
SiteSS 0.04306 0.02514 18.06310 1.713 0.104
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr)
SiteSS -0.293
eta_squared(Tol_Chl.lme_IN)
# Effect Size for ANOVA (Type III)
Parameter | Eta2 (partial) | 95% CI
-----------------------------------------
Site | 0.14 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_Chl_IN.res<-data.frame(summary(Tol_Chl.lme_IN)$coefficients)
names(Tol_Chl_IN.res)<-names(Chl_IN.res)[1:5]
Tol_Chl_IN.res$Predictor<-rep("Site", nrow(Tol_Chl_IN.res))
Tol_Chl_IN.res$Response<-rep("Chlorophyll Retention", nrow(Tol_Chl_IN.res))
Tol_Chl_IN.res$EtaSq<-c(eta_squared(Tol_Chl.lme_IN)$Eta2)
##Combine results
Chl_IN.res<-rbind(Chl_IN.res, Tol_Chl_IN.res[2,])
##Summary statistics by Site and Origin
IN_TolChl.sum<-summarySE(Therm_IN, measurevar="Chl.prop", groupvars=c("Site"), na.rm=TRUE)
##Plot Average Chlorophyll Retention across Treatments
IN_TolChl.plot<-ggplot(IN_TolChl.sum, aes(x=Site, y=Chl.prop, colour=Site)) +
scale_colour_manual(values=Site.colors.o)+
geom_errorbar(aes(ymin=Chl.prop-se, ymax=Chl.prop+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
geom_point(size=point.sz, position=position_dodge(width=0.5))+
theme_classic()+
theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"))+
labs(x="Site", y="Chlorophyll Retention", colour="Site")+
ylim(0, 0.35); IN_TolChl.plot
#Check normality
hist(Therm_IN$PAM.prop)
shapiro.test(Therm_IN$PAM.prop)
Shapiro-Wilk normality test
data: Therm_IN$PAM.prop
W = 0.8289, p-value = 0.001154
#Not normal
hist(Therm_IN$PAM.prop^2)
shapiro.test(Therm_IN$PAM.prop^2)
Shapiro-Wilk normality test
data: Therm_IN$PAM.prop^2
W = 0.86491, p-value = 0.005132
##Still not normal
##Model
#Function of Site, with Genotype as a Random effect
Tol_PAM.lme_IN<-lmer(PAM.prop~Site+(1|Genotype), data=Therm_IN)
##Check residuals
Tol_PAM.lme_res_IN <- simulateResiduals(fittedModel = Tol_PAM.lme_IN, plot = F)
plot(Tol_PAM.lme_res_IN)
##Model results
summary(Tol_PAM.lme_IN)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: PAM.prop ~ Site + (1 | Genotype)
Data: Therm_IN
REML criterion at convergence: -50.5
Scaled residuals:
Min 1Q Median 3Q Max
-2.7860 -0.2498 0.1225 0.5362 1.4591
Random effects:
Groups Name Variance Std.Dev.
Genotype (Intercept) 0.002221 0.04713
Residual 0.003538 0.05948
Number of obs: 23, groups: Genotype, 3
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.89138 0.03262 2.77551 27.325 0.000185 ***
SiteSS 0.05230 0.02487 19.00103 2.103 0.049053 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr)
SiteSS -0.399
eta_squared(Tol_PAM.lme_IN)
# Effect Size for ANOVA (Type III)
Parameter | Eta2 (partial) | 95% CI
-----------------------------------------
Site | 0.19 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_PAM_IN.res<-data.frame(summary(Tol_PAM.lme_IN)$coefficients)
names(Tol_PAM_IN.res)<-names(PAM_IN.res)[1:5]
Tol_PAM_IN.res$Predictor<-rep("Site", nrow(Tol_PAM_IN.res))
Tol_PAM_IN.res$Response<-rep("FvFm Retention", nrow(Tol_PAM_IN.res))
Tol_PAM_IN.res$EtaSq<-c(eta_squared(Tol_PAM.lme_IN)$Eta2)
##Combine results
PAM_IN.res<-rbind(PAM_IN.res, Tol_PAM_IN.res[2,])
Not normal, but model residuals are OK.
##Summary statistics by Site and Origin
IN_TolPAM.sum<-summarySE(Therm_IN, measurevar="PAM.prop", groupvars=c("Site"), na.rm=TRUE)
##Plot Average FvFm Retention across Treatments
IN_TolPAM.plot<-ggplot(IN_TolPAM.sum, aes(x=Site, y=PAM.prop, colour=Site)) +
scale_colour_manual(values=Site.colors.o)+
geom_errorbar(aes(ymin=PAM.prop-se, ymax=PAM.prop+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
geom_point(size=point.sz, position=position_dodge(width=0.5))+
theme_classic()+
theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"))+
labs(x="Site", y="Photochemical Efficiency Retention", colour="Site")+
ylim(0.6, 1)+
annotate("text", x=1.5, y=1, label="*", size=sig.sz, fontface="bold"); IN_TolPAM.plot
##Subset Timepoint 1
Therm_TP1<-subset(Therm_H, TimeP=="TP1")
#Check normality
hist(Therm_TP1$Sym.prop)
shapiro.test(Therm_TP1$Sym.prop)
Shapiro-Wilk normality test
data: Therm_TP1$Sym.prop
W = 0.91699, p-value = 0.002332
#Not normal
hist(log(Therm_TP1$Sym.prop+1))
shapiro.test(log(Therm_TP1$Sym.prop+1))
Shapiro-Wilk normality test
data: log(Therm_TP1$Sym.prop + 1)
W = 0.92258, p-value = 0.003671
##Still not normal
##Model with no transformation
#Function of Site and Origin, with Genotype as a Random effect
#Interaction between Site and Origin
Tol_Sym.lme_TP1<-lmer(Sym.prop~Origin*Site+(1|Genotype), data=Therm_TP1)
##Check residuals
Tol_Sym.lme_res_TP1 <- simulateResiduals(fittedModel = Tol_Sym.lme_TP1, plot = F)
plot(Tol_Sym.lme_res_TP1)
##Model results
summary(Tol_Sym.lme_TP1)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Sym.prop ~ Origin * Site + (1 | Genotype)
Data: Therm_TP1
REML criterion at convergence: -3.1
Scaled residuals:
Min 1Q Median 3Q Max
-3.2276 -0.5910 -0.0696 0.7547 1.8547
Random effects:
Groups Name Variance Std.Dev.
Genotype (Intercept) 0.07291 0.2700
Residual 0.03714 0.1927
Number of obs: 48, groups: Genotype, 3
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.35578 0.16552 2.38641 2.149 0.1436
OriginTransplant 0.14490 0.07867 42.00000 1.842 0.0726 .
SiteSS 0.06829 0.07867 42.00000 0.868 0.3903
OriginTransplant:SiteSS -0.15161 0.11126 42.00000 -1.363 0.1803
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) OrgnTr SiteSS
OrgnTrnspln -0.238
SiteSS -0.238 0.500
OrgnTrn:SSS 0.168 -0.707 -0.707
eta_squared(Tol_Sym.lme_TP1)
# Effect Size for ANOVA (Type III)
Parameter | Eta2 (partial) | 95% CI
-------------------------------------------
Origin | 0.04 | [0.00, 1.00]
Site | 4.34e-04 | [0.00, 1.00]
Origin:Site | 0.04 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_Sym_TP1.res<-data.frame(summary(Tol_Sym.lme_TP1)$coefficients[-1,])
Tol_Sym_TP1.res$Predictor<-c("Origin", "Site", "Origin x Site")
Tol_Sym_TP1.res$EtaSq<-c(eta_squared(Tol_Sym.lme_TP1)$Eta2)
Tol_Sym_TP1.res$Response<-rep("Symbiont Retention", nrow(Tol_Sym_TP1.res))
Not normal, but model residuals are OK.
##Pairwise Comparisons of Interest
#Native vs Transplant within a Site
Tol_Sym.pair_TP1<-emmeans(Tol_Sym.lme_TP1, pairwise~Origin | Site)
Tol_Sym.pair_TP1
$emmeans
Site = KL:
Origin emmean SE df lower.CL upper.CL
Native 0.356 0.166 2.39 -0.257 0.968
Transplant 0.501 0.166 2.39 -0.112 1.113
Site = SS:
Origin emmean SE df lower.CL upper.CL
Native 0.424 0.166 2.39 -0.188 1.036
Transplant 0.417 0.166 2.39 -0.195 1.030
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
Native - Transplant -0.14490 0.0787 42 -1.842 0.0726
Site = SS:
contrast estimate SE df t.ratio p.value
Native - Transplant 0.00671 0.0787 42 0.085 0.9325
Degrees-of-freedom method: kenward-roger
#Calculate standardized effect sizes for pairwise comparisons
Tol_Sym.pair.es_TP1<-data.frame(eff_size(Tol_Sym.pair_TP1, sigma=sigma(Tol_Sym.lme_TP1), edf=df.residual(Tol_Sym.lme_TP1)))
Since 'object' is a list, we are using the contrasts already present.
Tol_Sym.pair.es_TP1
Tol_Sym.pair.es_TP1<-Tol_Sym.pair.es_TP1 %>% dplyr::rename(contrast.se = contrast, SE.es = SE, df.es = df)
##Save Results
Tol_Sym.pair_TP1.res<-merge(data.frame(Tol_Sym.pair_TP1$contrasts), Tol_Sym.pair.es_TP1[,-c(1)])
Tol_Sym.pair_TP1.res$Response<-rep("Symbiont Retention", nrow(Tol_Sym.pair_TP1.res))
##Add Significance levels
Tol_Sym.pair_TP1.res$Sig<-ifelse(Tol_Sym.pair_TP1.res$p.value<0.001, "***",
ifelse(Tol_Sym.pair_TP1.res$p.value<0.01, "**",
ifelse(Tol_Sym.pair_TP1.res$p.value<0.05, "*",
ifelse(Tol_Sym.pair_TP1.res$p.value<0.1, "-", NA))))
Tol_Sym.pair_TP1.res
##Summary statistics by Site and Origin
TP1_TolSym.sum<-summarySE(Therm_TP1, measurevar="Sym.prop", groupvars=c("Site", "Origin", "Site.Orig"), na.rm=TRUE)
##Plot Average Symbiont Retention across Treatments
TP1_TolSym.plot<-ggplot(TP1_TolSym.sum, aes(x=Site, y=Sym.prop, colour=Site.Orig)) +
scale_colour_manual(values=Orig.colors.o)+
geom_errorbar(aes(ymin=Sym.prop-se, ymax=Sym.prop+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
geom_point(size=point.sz, position=position_dodge(width=0.5))+
theme_classic()+
theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"))+
labs(x="Site and Origin", y="Symbiont Retention", colour="Origin")+
ylim(0, 1)+
annotate("text", x=1, y=0.65, label="-", size=sig.sz, fontface="bold"); TP1_TolSym.plot
##Pairwise Comparisons of Interest
#Native vs Transplant within a Site and Timepoint
Tol_Chl.pair<-emmeans(Tol_Chl.lme, pairwise~Origin | Site*TimeP)
Tol_Chl.pair
$emmeans
Site = KL, TimeP = TP1:
Origin emmean SE df lower.CL upper.CL
Native 0.1556 0.0408 41.3 0.07319 0.238
Transplant 0.1697 0.0408 41.3 0.08727 0.252
Site = SS, TimeP = TP1:
Origin emmean SE df lower.CL upper.CL
Native 0.1849 0.0408 41.3 0.10249 0.267
Transplant 0.1726 0.0408 41.3 0.09010 0.255
Site = KL, TimeP = TP2:
Origin emmean SE df lower.CL upper.CL
Native 0.3108 0.0408 41.3 0.22832 0.393
Transplant 0.2723 0.0408 41.3 0.18988 0.355
Site = SS, TimeP = TP2:
Origin emmean SE df lower.CL upper.CL
Native 0.2456 0.0408 41.3 0.16318 0.328
Transplant 0.2688 0.0408 41.3 0.18639 0.351
Site = KL, TimeP = TP3:
Origin emmean SE df lower.CL upper.CL
Native 0.5549 0.0408 41.3 0.47243 0.637
Transplant 0.6116 0.0408 41.3 0.52910 0.694
Site = SS, TimeP = TP3:
Origin emmean SE df lower.CL upper.CL
Native 0.6275 0.0408 41.3 0.54508 0.710
Transplant 0.6382 0.0408 41.3 0.55573 0.721
Site = KL, TimeP = TP4:
Origin emmean SE df lower.CL upper.CL
Native 0.1000 0.0408 41.3 0.01759 0.182
Transplant 0.0858 0.0408 41.3 0.00330 0.168
Site = SS, TimeP = TP4:
Origin emmean SE df lower.CL upper.CL
Native 0.0723 0.0424 46.5 -0.01305 0.158
Transplant 0.0895 0.0408 41.3 0.00707 0.172
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
$contrasts
Site = KL, TimeP = TP1:
contrast estimate SE df t.ratio p.value
Native - Transplant -0.0141 0.0533 173 -0.264 0.7919
Site = SS, TimeP = TP1:
contrast estimate SE df t.ratio p.value
Native - Transplant 0.0124 0.0533 173 0.233 0.8163
Site = KL, TimeP = TP2:
contrast estimate SE df t.ratio p.value
Native - Transplant 0.0384 0.0533 173 0.722 0.4715
Site = SS, TimeP = TP2:
contrast estimate SE df t.ratio p.value
Native - Transplant -0.0232 0.0533 173 -0.436 0.6635
Site = KL, TimeP = TP3:
contrast estimate SE df t.ratio p.value
Native - Transplant -0.0567 0.0533 173 -1.064 0.2889
Site = SS, TimeP = TP3:
contrast estimate SE df t.ratio p.value
Native - Transplant -0.0107 0.0533 173 -0.200 0.8417
Site = KL, TimeP = TP4:
contrast estimate SE df t.ratio p.value
Native - Transplant 0.0143 0.0533 173 0.268 0.7888
Site = SS, TimeP = TP4:
contrast estimate SE df t.ratio p.value
Native - Transplant -0.0172 0.0545 173 -0.316 0.7523
Degrees-of-freedom method: kenward-roger
#Calculate standardized effect sizes for pairwise comparisons
Tol_Chl.pair.es<-data.frame(eff_size(Tol_Chl.pair, sigma=sigma(Tol_Chl.lme), edf=df.residual(Tol_Chl.lme)))
Since 'object' is a list, we are using the contrasts already present.
Tol_Chl.pair.es
Tol_Chl.pair.es<-Tol_Chl.pair.es %>% dplyr::rename(contrast.se = contrast, SE.es = SE, df.es = df)
##Save Results
Tol_Chl.pair.res<-merge(data.frame(Tol_Chl.pair$contrasts), Tol_Chl.pair.es[,-c(1)])
Tol_Chl.pair.res$Response<-rep("Chlorophyll Retention", nrow(Tol_Chl.pair.res))
Not normal, but model residuals are OK.
##Pairwise Comparisons of Interest
#Native vs Transplant within a Site
Tol_Chl.pair_TP1<-emmeans(Tol_Chl.lme_TP1, pairwise~Origin | Site)
Tol_Chl.pair_TP1
$emmeans
Site = KL:
Origin emmean SE df lower.CL upper.CL
Native 0.156 0.0657 2.27 -0.0970 0.408
Transplant 0.170 0.0657 2.27 -0.0829 0.422
Site = SS:
Origin emmean SE df lower.CL upper.CL
Native 0.185 0.0657 2.27 -0.0677 0.438
Transplant 0.173 0.0657 2.27 -0.0801 0.425
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
Native - Transplant -0.0141 0.0266 42 -0.529 0.5997
Site = SS:
contrast estimate SE df t.ratio p.value
Native - Transplant 0.0124 0.0266 42 0.466 0.6439
Degrees-of-freedom method: kenward-roger
#Calculate standardized effect sizes for pairwise comparisons
Tol_Chl.pair.es_TP1<-data.frame(eff_size(Tol_Chl.pair_TP1, sigma=sigma(Tol_Chl.lme_TP1), edf=df.residual(Tol_Chl.lme_TP1)))
Since 'object' is a list, we are using the contrasts already present.
Tol_Chl.pair.es_TP1
Tol_Chl.pair.es_TP1<-Tol_Chl.pair.es_TP1 %>% dplyr::rename(contrast.se = contrast, SE.es = SE, df.es = df)
##Save Results
Tol_Chl.pair_TP1.res<-merge(data.frame(Tol_Chl.pair_TP1$contrasts), Tol_Chl.pair.es_TP1[,-c(1)])
Tol_Chl.pair_TP1.res$Response<-rep("Chlorophyll Retention", nrow(Tol_Chl.pair_TP1.res))
##Add Significance levels
Tol_Chl.pair_TP1.res$Sig<-ifelse(Tol_Chl.pair_TP1.res$p.value<0.001, "***",
ifelse(Tol_Chl.pair_TP1.res$p.value<0.01, "**",
ifelse(Tol_Chl.pair_TP1.res$p.value<0.05, "*",
ifelse(Tol_Chl.pair_TP1.res$p.value<0.1, "-", NA))))
Tol_Chl.pair_TP1.res
##Summary statistics by Site and Origin
TP1_TolChl.sum<-summarySE(Therm_TP1, measurevar="Chl.prop", groupvars=c("Site", "Origin", "Site.Orig"), na.rm=TRUE)
##Plot Average Chlorophyll Retention across Treatments
TP1_TolChl.plot<-ggplot(TP1_TolChl.sum, aes(x=Site, y=Chl.prop, colour=Site.Orig)) +
scale_colour_manual(values=Orig.colors.o)+
geom_errorbar(aes(ymin=Chl.prop-se, ymax=Chl.prop+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
geom_point(size=point.sz, position=position_dodge(width=0.5))+
theme_classic()+
theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"))+
labs(x="Site and Origin", y="Chlorophyll Retention", colour="Origin")+
ylim(0, 0.35); TP1_TolChl.plot
#Check normality
hist(Therm_TP1$PAM.prop)
shapiro.test(Therm_TP1$PAM.prop)
Shapiro-Wilk normality test
data: Therm_TP1$PAM.prop
W = 0.84073, p-value = 1.245e-05
#Not normal
hist(Therm_TP1$PAM.prop^2)
shapiro.test(Therm_TP1$PAM.prop^2)
Shapiro-Wilk normality test
data: Therm_TP1$PAM.prop^2
W = 0.88155, p-value = 0.0001681
##Still not normal
##Model with no transformation
#Function of Site and Origin, with Genotype as a Random effect
#Interaction between Site and Origin
Tol_PAM.lme_TP1<-lmer(PAM.prop~Origin*Site+(1|Genotype), data=Therm_TP1)
##Check residuals
Tol_PAM.lme_res_TP1 <- simulateResiduals(fittedModel = Tol_PAM.lme_TP1, plot = F)
plot(Tol_PAM.lme_res_TP1)
##Model results
summary(Tol_PAM.lme_TP1)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: PAM.prop ~ Origin * Site + (1 | Genotype)
Data: Therm_TP1
REML criterion at convergence: -84.6
Scaled residuals:
Min 1Q Median 3Q Max
-2.3058 -0.3498 0.1733 0.6717 1.7818
Random effects:
Groups Name Variance Std.Dev.
Genotype (Intercept) 0.002735 0.05230
Residual 0.006214 0.07883
Number of obs: 48, groups: Genotype, 3
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.86181 0.03781 3.74555 22.794 3.72e-05 ***
OriginTransplant 0.06813 0.03218 42.00000 2.117 0.0402 *
SiteSS -0.01277 0.03218 42.00000 -0.397 0.6936
OriginTransplant:SiteSS -0.07942 0.04551 42.00000 -1.745 0.0883 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) OrgnTr SiteSS
OrgnTrnspln -0.426
SiteSS -0.426 0.500
OrgnTrn:SSS 0.301 -0.707 -0.707
eta_squared(Tol_PAM.lme_TP1)
# Effect Size for ANOVA (Type III)
Parameter | Eta2 (partial) | 95% CI
-------------------------------------------
Origin | 0.04 | [0.00, 1.00]
Site | 0.11 | [0.01, 1.00]
Origin:Site | 0.07 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_PAM_TP1.res<-data.frame(summary(Tol_PAM.lme_TP1)$coefficients[-1,])
Tol_PAM_TP1.res$Predictor<-c("Origin", "Site", "Origin x Site")
Tol_PAM_TP1.res$EtaSq<-c(eta_squared(Tol_PAM.lme_TP1)$Eta2)
Tol_PAM_TP1.res$Response<-rep("FvFm Retention", nrow(Tol_PAM_TP1.res))
Not normal, but model residuals are OK.
##Pairwise Comparisons of Interest
#Native vs Transplant within a Site
Tol_PAM.pair_TP1<-emmeans(Tol_PAM.lme_TP1, pairwise~Origin | Site)
Tol_PAM.pair_TP1
$emmeans
Site = KL:
Origin emmean SE df lower.CL upper.CL
Native 0.862 0.0378 3.75 0.754 0.970
Transplant 0.930 0.0378 3.75 0.822 1.038
Site = SS:
Origin emmean SE df lower.CL upper.CL
Native 0.849 0.0378 3.75 0.741 0.957
Transplant 0.838 0.0378 3.75 0.730 0.946
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
Native - Transplant -0.0681 0.0322 42 -2.117 0.0402
Site = SS:
contrast estimate SE df t.ratio p.value
Native - Transplant 0.0113 0.0322 42 0.351 0.7276
Degrees-of-freedom method: kenward-roger
#Calculate standardized effect sizes for pairwise comparisons
Tol_PAM.pair.es_TP1<-data.frame(eff_size(Tol_PAM.pair_TP1, sigma=sigma(Tol_PAM.lme_TP1), edf=df.residual(Tol_PAM.lme_TP1)))
Since 'object' is a list, we are using the contrasts already present.
Tol_PAM.pair.es_TP1
Tol_PAM.pair.es_TP1<-Tol_PAM.pair.es_TP1 %>% dplyr::rename(contrast.se = contrast, SE.es = SE, df.es = df)
##Save Results
Tol_PAM.pair_TP1.res<-merge(data.frame(Tol_PAM.pair_TP1$contrasts), Tol_PAM.pair.es_TP1[,-c(1)])
Tol_PAM.pair_TP1.res$Response<-rep("FvFm Retention", nrow(Tol_PAM.pair_TP1.res))
##Add Significance levels
Tol_PAM.pair_TP1.res$Sig<-ifelse(Tol_PAM.pair_TP1.res$p.value<0.001, "***",
ifelse(Tol_PAM.pair_TP1.res$p.value<0.01, "**",
ifelse(Tol_PAM.pair_TP1.res$p.value<0.05, "*",
ifelse(Tol_PAM.pair_TP1.res$p.value<0.1, "-", NA))))
Tol_PAM.pair_TP1.res
##Summary statistics by Site and Origin
TP1_TolPAM.sum<-summarySE(Therm_TP1, measurevar="PAM.prop", groupvars=c("Site", "Origin", "Site.Orig"), na.rm=TRUE)
##Plot Average FvFm Retention across Treatments
TP1_TolPAM.plot<-ggplot(TP1_TolPAM.sum, aes(x=Site, y=PAM.prop, colour=Site.Orig)) +
scale_colour_manual(values=Orig.colors.o)+
geom_errorbar(aes(ymin=PAM.prop-se, ymax=PAM.prop+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
geom_point(size=point.sz, position=position_dodge(width=0.5))+
theme_classic()+
theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"))+
labs(x="Site and Origin", y="Photochemical Efficiency Retention", colour="Origin")+
ylim(0.6, 1)+
annotate("text", x=1, y=.96, label="*", size=sig.sz, fontface="bold"); TP1_TolPAM.plot
##Combine Results
Tol_TP1.res<-rbind(Tol_PAM_TP1.res, Tol_Sym_TP1.res, Tol_Chl_TP1.res)
Tol.pair.res_TP1<-rbind(Tol_PAM.pair_TP1.res, Tol_Sym.pair_TP1.res, Tol_Chl.pair_TP1.res)
##Add Timepoint
Tol_TP1.res$TimeP<-rep("TP1", nrow(Tol_TP1.res))
Tol.pair.res_TP1$TimeP<-rep("TP1", nrow(Tol.pair.res_TP1))
##Subset Timepoint 2
Therm_TP2<-subset(Therm_H, TimeP=="TP2")
#Check normality
hist(Therm_TP2$Sym.prop)
shapiro.test(Therm_TP2$Sym.prop)
Shapiro-Wilk normality test
data: Therm_TP2$Sym.prop
W = 0.89527, p-value = 0.0004439
#Not normal
hist(Therm_TP2$Sym.prop^2)
shapiro.test(Therm_TP2$Sym.prop^2)
Shapiro-Wilk normality test
data: Therm_TP2$Sym.prop^2
W = 0.8757, p-value = 0.0001129
##Still not normal
##Model with no transformation
#Function of Site and Origin, with Genotype as a Random effect
#Interaction between Site and Origin
Tol_Sym.lme_TP2<-lmer(Sym.prop~Origin*Site+(1|Genotype), data=Therm_TP2)
##Check residuals
Tol_Sym.lme_res_TP2 <- simulateResiduals(fittedModel = Tol_Sym.lme_TP2, plot = F)
plot(Tol_Sym.lme_res_TP2)
##Model results
summary(Tol_Sym.lme_TP2)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Sym.prop ~ Origin * Site + (1 | Genotype)
Data: Therm_TP2
REML criterion at convergence: -3.4
Scaled residuals:
Min 1Q Median 3Q Max
-2.20381 -0.54227 0.01202 0.47665 2.09092
Random effects:
Groups Name Variance Std.Dev.
Genotype (Intercept) 0.04681 0.2164
Residual 0.03767 0.1941
Number of obs: 48, groups: Genotype, 3
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.805833 0.136905 2.613359 5.886 0.0142 *
OriginTransplant -0.107900 0.079237 42.000000 -1.362 0.1805
SiteSS -0.178425 0.079237 42.000000 -2.252 0.0296 *
OriginTransplant:SiteSS 0.005183 0.112058 42.000000 0.046 0.9633
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) OrgnTr SiteSS
OrgnTrnspln -0.289
SiteSS -0.289 0.500
OrgnTrn:SSS 0.205 -0.707 -0.707
eta_squared(Tol_Sym.lme_TP2)
# Effect Size for ANOVA (Type III)
Parameter | Eta2 (partial) | 95% CI
-------------------------------------------
Origin | 0.08 | [0.00, 1.00]
Site | 0.19 | [0.04, 1.00]
Origin:Site | 5.09e-05 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_Sym_TP2.res<-data.frame(summary(Tol_Sym.lme_TP2)$coefficients[-1,])
Tol_Sym_TP2.res$Predictor<-c("Origin", "Site", "Origin x Site")
Tol_Sym_TP2.res$EtaSq<-c(eta_squared(Tol_Sym.lme_TP2)$Eta2)
Tol_Sym_TP2.res$Response<-rep("Symbiont Retention", nrow(Tol_Sym_TP2.res))
Not normal, but model residuals are OK.
##Pairwise Comparisons of Interest
#Native vs Transplant within a Site
Tol_Sym.pair_TP2<-emmeans(Tol_Sym.lme_TP2, pairwise~Origin | Site)
Tol_Sym.pair_TP2
$emmeans
Site = KL:
Origin emmean SE df lower.CL upper.CL
Native 0.806 0.137 2.61 0.3313 1.280
Transplant 0.698 0.137 2.61 0.2234 1.172
Site = SS:
Origin emmean SE df lower.CL upper.CL
Native 0.627 0.137 2.61 0.1529 1.102
Transplant 0.525 0.137 2.61 0.0502 0.999
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
Native - Transplant 0.108 0.0792 42 1.362 0.1805
Site = SS:
contrast estimate SE df t.ratio p.value
Native - Transplant 0.103 0.0792 42 1.296 0.2019
Degrees-of-freedom method: kenward-roger
#Calculate standardized effect sizes for pairwise comparisons
Tol_Sym.pair.es_TP2<-data.frame(eff_size(Tol_Sym.pair_TP2, sigma=sigma(Tol_Sym.lme_TP2), edf=df.residual(Tol_Sym.lme_TP2)))
Since 'object' is a list, we are using the contrasts already present.
Tol_Sym.pair.es_TP2
Tol_Sym.pair.es_TP2<-Tol_Sym.pair.es_TP2 %>% dplyr::rename(contrast.se = contrast, SE.es = SE, df.es = df)
##Save Results
Tol_Sym.pair_TP2.res<-merge(data.frame(Tol_Sym.pair_TP2$contrasts), Tol_Sym.pair.es_TP2[,-c(1)])
Tol_Sym.pair_TP2.res$Response<-rep("Symbiont Retention", nrow(Tol_Sym.pair_TP2.res))
##Add Significance levels
Tol_Sym.pair_TP2.res$Sig<-ifelse(Tol_Sym.pair_TP2.res$p.value<0.001, "***",
ifelse(Tol_Sym.pair_TP2.res$p.value<0.01, "**",
ifelse(Tol_Sym.pair_TP2.res$p.value<0.05, "*",
ifelse(Tol_Sym.pair_TP2.res$p.value<0.1, "-", NA))))
Tol_Sym.pair_TP2.res
##Summary statistics by Site and Origin
TP2_TolSym.sum<-summarySE(Therm_TP2, measurevar="Sym.prop", groupvars=c("Site", "Origin", "Site.Orig"), na.rm=TRUE)
##Plot Average Retention across Treatments
TP2_TolSym.plot<-ggplot(TP2_TolSym.sum, aes(x=Site, y=Sym.prop, colour=Site.Orig)) +
scale_colour_manual(values=Orig.colors.o)+
geom_errorbar(aes(ymin=Sym.prop-se, ymax=Sym.prop+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
geom_point(size=point.sz, position=position_dodge(width=0.5))+
theme_classic()+
theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"))+
labs(x="Site and Origin", y="Symbiont Retention", colour="Origin")+
ylim(0, 1)+
geom_bracket(xmin=1, xmax=2, y.position=.04, label.size=levels.sz, label="*", inherit.aes = FALSE); TP2_TolSym.plot
#Check normality
hist(Therm_TP2$Chl.prop)
shapiro.test(Therm_TP2$Chl.prop)
Shapiro-Wilk normality test
data: Therm_TP2$Chl.prop
W = 0.97677, p-value = 0.4523
#Normal
##Model
#Function of Site and Origin, with Genotype as a Random effect
#Interaction between Site and Origin
Tol_Chl.lme_TP2<-lmer(Chl.prop~Origin*Site+(1|Genotype), data=Therm_TP2)
##Check residuals
Tol_Chl.lme_res_TP2 <- simulateResiduals(fittedModel = Tol_Chl.lme_TP2, plot = F)
plot(Tol_Chl.lme_res_TP2)
##Model results
summary(Tol_Chl.lme_TP2)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: Chl.prop ~ Origin * Site + (1 | Genotype)
Data: Therm_TP2
REML criterion at convergence: -81
Scaled residuals:
Min 1Q Median 3Q Max
-2.0024 -0.7622 -0.1174 0.6824 2.9666
Random effects:
Groups Name Variance Std.Dev.
Genotype (Intercept) 0.006146 0.07840
Residual 0.006538 0.08086
Number of obs: 48, groups: Genotype, 3
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.31077 0.05093 2.81331 6.102 0.0106 *
OriginTransplant -0.03844 0.03301 42.00000 -1.165 0.2508
SiteSS -0.06514 0.03301 42.00000 -1.973 0.0550 .
OriginTransplant:SiteSS 0.06166 0.04668 42.00000 1.321 0.1937
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) OrgnTr SiteSS
OrgnTrnspln -0.324
SiteSS -0.324 0.500
OrgnTrn:SSS 0.229 -0.707 -0.707
eta_squared(Tol_Chl.lme_TP2)
# Effect Size for ANOVA (Type III)
Parameter | Eta2 (partial) | 95% CI
-------------------------------------------
Origin | 2.53e-03 | [0.00, 1.00]
Site | 0.05 | [0.00, 1.00]
Origin:Site | 0.04 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_Chl_TP2.res<-data.frame(summary(Tol_Chl.lme_TP2)$coefficients[-1,])
Tol_Chl_TP2.res$Predictor<-c("Origin", "Site", "Origin x Site")
Tol_Chl_TP2.res$EtaSq<-c(eta_squared(Tol_Chl.lme_TP2)$Eta2)
Tol_Chl_TP2.res$Response<-rep("Chlorophyll Retention", nrow(Tol_Chl_TP2.res))
##Pairwise Comparisons of Interest
#Native vs Transplant within a Site
Tol_Chl.pair_TP2<-emmeans(Tol_Chl.lme_TP2, pairwise~Origin | Site)
Tol_Chl.pair_TP2
$emmeans
Site = KL:
Origin emmean SE df lower.CL upper.CL
Native 0.311 0.0509 2.81 0.1424 0.479
Transplant 0.272 0.0509 2.81 0.1040 0.441
Site = SS:
Origin emmean SE df lower.CL upper.CL
Native 0.246 0.0509 2.81 0.0773 0.414
Transplant 0.269 0.0509 2.81 0.1005 0.437
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
Native - Transplant 0.0384 0.033 42 1.165 0.2508
Site = SS:
contrast estimate SE df t.ratio p.value
Native - Transplant -0.0232 0.033 42 -0.703 0.4857
Degrees-of-freedom method: kenward-roger
#Calculate standardized effect sizes for pairwise comparisons
Tol_Chl.pair.es_TP2<-data.frame(eff_size(Tol_Chl.pair_TP2, sigma=sigma(Tol_Chl.lme_TP2), edf=df.residual(Tol_Chl.lme_TP2)))
Since 'object' is a list, we are using the contrasts already present.
Tol_Chl.pair.es_TP2
Tol_Chl.pair.es_TP2<-Tol_Chl.pair.es_TP2 %>% dplyr::rename(contrast.se = contrast, SE.es = SE, df.es = df)
##Save Results
Tol_Chl.pair_TP2.res<-merge(data.frame(Tol_Chl.pair_TP2$contrasts), Tol_Chl.pair.es_TP2[,-c(1)])
Tol_Chl.pair_TP2.res$Response<-rep("Chlorophyll Retention", nrow(Tol_Chl.pair_TP2.res))
##Add Significance levels
Tol_Chl.pair_TP2.res$Sig<-ifelse(Tol_Chl.pair_TP2.res$p.value<0.001, "***",
ifelse(Tol_Chl.pair_TP2.res$p.value<0.01, "**",
ifelse(Tol_Chl.pair_TP2.res$p.value<0.05, "*",
ifelse(Tol_Chl.pair_TP2.res$p.value<0.1, "-", NA))))
Tol_Chl.pair_TP2.res
##Summary statistics by Site and Origin
TP2_TolChl.sum<-summarySE(Therm_TP2, measurevar="Chl.prop", groupvars=c("Site", "Origin", "Site.Orig"), na.rm=TRUE)
##Plot Average Retention across Treatments
TP2_TolChl.plot<-ggplot(TP2_TolChl.sum, aes(x=Site, y=Chl.prop, colour=Site.Orig)) +
scale_colour_manual(values=Orig.colors.o)+
geom_errorbar(aes(ymin=Chl.prop-se, ymax=Chl.prop+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
geom_point(size=point.sz, position=position_dodge(width=0.5))+
theme_classic()+
theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"))+
labs(x="Site and Origin", y="Chlorophyll Retention", colour="Origin")+
ylim(0, 0.35)+
geom_bracket(xmin=1, xmax=2, y.position=.02, label.size=levels.sz, label="-", inherit.aes = FALSE); TP2_TolChl.plot
#Check normality
hist(Therm_TP2$PAM.prop)
shapiro.test(Therm_TP2$PAM.prop)
Shapiro-Wilk normality test
data: Therm_TP2$PAM.prop
W = 0.73527, p-value = 5.896e-08
#Not normal
hist(Therm_TP2$PAM.prop^2)
shapiro.test(Therm_TP2$PAM.prop^2)
Shapiro-Wilk normality test
data: Therm_TP2$PAM.prop^2
W = 0.80435, p-value = 1.633e-06
##Still not normal
##Model
#Function of Site and Origin, with Genotype as a Random effect
#Interaction between Site and Origin
Tol_PAM.lme_TP2<-lmer(PAM.prop~Origin*Site+(1|Genotype), data=Therm_TP2)
##Check residuals
Tol_PAM.lme_res_TP2 <- simulateResiduals(fittedModel = Tol_PAM.lme_TP2, plot = F)
plot(Tol_PAM.lme_res_TP2)
##Model results
summary(Tol_PAM.lme_TP2)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: PAM.prop ~ Origin * Site + (1 | Genotype)
Data: Therm_TP2
REML criterion at convergence: -102.1
Scaled residuals:
Min 1Q Median 3Q Max
-3.8326 -0.3753 -0.0042 0.5350 1.7496
Random effects:
Groups Name Variance Std.Dev.
Genotype (Intercept) 0.001675 0.04093
Residual 0.004192 0.06475
Number of obs: 48, groups: Genotype, 3
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.85622 0.03013 3.92146 28.419 1.09e-05 ***
OriginTransplant 0.07526 0.02643 42.00000 2.847 0.00680 **
SiteSS 0.06487 0.02643 42.00000 2.454 0.01834 *
OriginTransplant:SiteSS -0.10206 0.03738 42.00000 -2.730 0.00921 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) OrgnTr SiteSS
OrgnTrnspln -0.439
SiteSS -0.439 0.500
OrgnTrn:SSS 0.310 -0.707 -0.707
eta_squared(Tol_PAM.lme_TP2)
# Effect Size for ANOVA (Type III)
Parameter | Eta2 (partial) | 95% CI
-------------------------------------------
Origin | 0.04 | [0.00, 1.00]
Site | 0.01 | [0.00, 1.00]
Origin:Site | 0.15 | [0.02, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_PAM_TP2.res<-data.frame(summary(Tol_PAM.lme_TP2)$coefficients[-1,])
Tol_PAM_TP2.res$Predictor<-c("Origin", "Site", "Origin x Site")
Tol_PAM_TP2.res$EtaSq<-c(eta_squared(Tol_PAM.lme_TP2)$Eta2)
Tol_PAM_TP2.res$Response<-rep("FvFm Retention", nrow(Tol_PAM_TP2.res))
Not normal, but model residuals are OK.
##Pairwise Comparisons of Interest
#Native vs Transplant within a Site
Tol_PAM.pair_TP2<-emmeans(Tol_PAM.lme_TP2, pairwise~Origin | Site)
Tol_PAM.pair_TP2
$emmeans
Site = KL:
Origin emmean SE df lower.CL upper.CL
Native 0.856 0.0301 3.92 0.772 0.941
Transplant 0.931 0.0301 3.92 0.847 1.016
Site = SS:
Origin emmean SE df lower.CL upper.CL
Native 0.921 0.0301 3.92 0.837 1.005
Transplant 0.894 0.0301 3.92 0.810 0.979
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
Native - Transplant -0.0753 0.0264 42 -2.847 0.0068
Site = SS:
contrast estimate SE df t.ratio p.value
Native - Transplant 0.0268 0.0264 42 1.014 0.3164
Degrees-of-freedom method: kenward-roger
#Calculate standardized effect sizes for pairwise comparisons
Tol_PAM.pair.es_TP2<-data.frame(eff_size(Tol_PAM.pair_TP2, sigma=sigma(Tol_PAM.lme_TP2), edf=df.residual(Tol_PAM.lme_TP2)))
Since 'object' is a list, we are using the contrasts already present.
Tol_PAM.pair.es_TP2
Tol_PAM.pair.es_TP2<-Tol_PAM.pair.es_TP2 %>% dplyr::rename(contrast.se = contrast, SE.es = SE, df.es = df)
##Save Results
Tol_PAM.pair_TP2.res<-merge(data.frame(Tol_PAM.pair_TP2$contrasts), Tol_PAM.pair.es_TP2[,-c(1)])
Tol_PAM.pair_TP2.res$Response<-rep("FvFm Retention", nrow(Tol_PAM.pair_TP2.res))
##Add Significance levels
Tol_PAM.pair_TP2.res$Sig<-ifelse(Tol_PAM.pair_TP2.res$p.value<0.001, "***",
ifelse(Tol_PAM.pair_TP2.res$p.value<0.01, "**",
ifelse(Tol_PAM.pair_TP2.res$p.value<0.05, "*",
ifelse(Tol_PAM.pair_TP2.res$p.value<0.1, "-", NA))))
Tol_PAM.pair_TP2.res
##Combine Results
Tol_TP2.res<-rbind(Tol_PAM_TP2.res, Tol_Sym_TP2.res, Tol_Chl_TP2.res)
Tol.pair.res_TP2<-rbind(Tol_PAM.pair_TP2.res, Tol_Sym.pair_TP2.res, Tol_Chl.pair_TP2.res)
##Add Timepoint
Tol_TP2.res$TimeP<-rep("TP2", nrow(Tol_TP2.res))
Tol.pair.res_TP2$TimeP<-rep("TP2", nrow(Tol.pair.res_TP2))
##Subset Timepoint 3
Therm_TP3<-subset(Therm_H, TimeP=="TP3")
#Check normality
hist(Therm_TP3$Sym.prop)
shapiro.test(Therm_TP3$Sym.prop)
Shapiro-Wilk normality test
data: Therm_TP3$Sym.prop
W = 0.78653, p-value = 6.531e-07
#Not normal
hist(Therm_TP3$Sym.prop^2)
shapiro.test(Therm_TP3$Sym.prop^2)
Shapiro-Wilk normality test
data: Therm_TP3$Sym.prop^2
W = 0.78554, p-value = 6.213e-07
##Still not normal but less skewed
##Model with no transformation
#Function of Site and Origin, with Genotype as a Random effect
#Interaction between Site and Origin
Tol_Sym.lme_TP3<-lmer(Sym.prop~Origin*Site+(1|Genotype), data=Therm_TP3)
##Check residuals
Tol_Sym.lme_res_TP3 <- simulateResiduals(fittedModel = Tol_Sym.lme_TP3, plot = F)
plot(Tol_Sym.lme_res_TP3)
##Model results
summary(Tol_Sym.lme_TP3)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Sym.prop ~ Origin * Site + (1 | Genotype)
Data: Therm_TP3
REML criterion at convergence: -54.7
Scaled residuals:
Min 1Q Median 3Q Max
-2.1418 -0.5595 0.3785 0.6079 1.6260
Random effects:
Groups Name Variance Std.Dev.
Genotype (Intercept) 0.0003012 0.01735
Residual 0.0132967 0.11531
Number of obs: 48, groups: Genotype, 3
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.82201 0.03476 16.65879 23.646 2.98e-14 ***
OriginTransplant 0.07462 0.04708 42.00000 1.585 0.1204
SiteSS 0.10814 0.04708 42.00000 2.297 0.0267 *
OriginTransplant:SiteSS -0.05670 0.06658 42.00000 -0.852 0.3992
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) OrgnTr SiteSS
OrgnTrnspln -0.677
SiteSS -0.677 0.500
OrgnTrn:SSS 0.479 -0.707 -0.707
eta_squared(Tol_Sym.lme_TP3)
# Effect Size for ANOVA (Type III)
Parameter | Eta2 (partial) | 95% CI
-------------------------------------------
Origin | 0.04 | [0.00, 1.00]
Site | 0.12 | [0.01, 1.00]
Origin:Site | 0.02 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_Sym_TP3.res<-data.frame(summary(Tol_Sym.lme_TP3)$coefficients[-1,])
Tol_Sym_TP3.res$Predictor<-c("Origin", "Site", "Origin x Site")
Tol_Sym_TP3.res$EtaSq<-c(eta_squared(Tol_Sym.lme_TP3)$Eta2)
Tol_Sym_TP3.res$Response<-rep("Symbiont Retention", nrow(Tol_Sym_TP3.res))
Not normal, but model residuals are OK.
##Pairwise Comparisons of Interest
#Native vs Transplant within a Site
Tol_Sym.pair_TP3<-emmeans(Tol_Sym.lme_TP3, pairwise~Origin | Site)
Tol_Sym.pair_TP3
$emmeans
Site = KL:
Origin emmean SE df lower.CL upper.CL
Native 0.822 0.0348 16.7 0.749 0.895
Transplant 0.897 0.0348 16.7 0.823 0.970
Site = SS:
Origin emmean SE df lower.CL upper.CL
Native 0.930 0.0348 16.7 0.857 1.004
Transplant 0.948 0.0348 16.7 0.875 1.022
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
Native - Transplant -0.0746 0.0471 42 -1.585 0.1204
Site = SS:
contrast estimate SE df t.ratio p.value
Native - Transplant -0.0179 0.0471 42 -0.381 0.7053
Degrees-of-freedom method: kenward-roger
#Calculate standardized effect sizes for pairwise comparisons
Tol_Sym.pair.es_TP3<-data.frame(eff_size(Tol_Sym.pair_TP3, sigma=sigma(Tol_Sym.lme_TP3), edf=df.residual(Tol_Sym.lme_TP3)))
Since 'object' is a list, we are using the contrasts already present.
Tol_Sym.pair.es_TP3
Tol_Sym.pair.es_TP3<-Tol_Sym.pair.es_TP3 %>% dplyr::rename(contrast.se = contrast, SE.es = SE, df.es = df)
##Save Results
Tol_Sym.pair_TP3.res<-merge(data.frame(Tol_Sym.pair_TP3$contrasts), Tol_Sym.pair.es_TP3[,-c(1)])
Tol_Sym.pair_TP3.res$Response<-rep("Symbiont Retention", nrow(Tol_Sym.pair_TP3.res))
##Add Significance levels
Tol_Sym.pair_TP3.res$Sig<-ifelse(Tol_Sym.pair_TP3.res$p.value<0.001, "***",
ifelse(Tol_Sym.pair_TP3.res$p.value<0.01, "**",
ifelse(Tol_Sym.pair_TP3.res$p.value<0.05, "*",
ifelse(Tol_Sym.pair_TP3.res$p.value<0.1, "-", NA))))
Tol_Sym.pair_TP3.res
##Summary statistics by Site and Origin
TP3_TolSym.sum<-summarySE(Therm_TP3, measurevar="Sym.prop", groupvars=c("Site", "Origin", "Site.Orig"), na.rm=TRUE)
##Plot Average Retention across Treatments
TP3_TolSym.plot<-ggplot(TP3_TolSym.sum, aes(x=Site, y=Sym.prop, colour=Site.Orig)) +
scale_colour_manual(values=Orig.colors.o)+
geom_errorbar(aes(ymin=Sym.prop-se, ymax=Sym.prop+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
geom_point(size=point.sz, position=position_dodge(width=0.5))+
theme_classic()+
theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"))+
labs(x="Site and Origin", y="Symbiont Retention", colour="Origin")+
ylim(0, 1)+
geom_bracket(xmin=1, xmax=2, y.position=.04, label.size=levels.sz, label="*", inherit.aes = FALSE); TP3_TolSym.plot
#Check normality
hist(Therm_TP3$Chl.prop)
shapiro.test(Therm_TP3$Chl.prop)
Shapiro-Wilk normality test
data: Therm_TP3$Chl.prop
W = 0.89534, p-value = 0.0004461
#Not normal
hist(log(Therm_TP3$Chl.prop+1))
shapiro.test(log(Therm_TP3$Chl.prop+1))
Shapiro-Wilk normality test
data: log(Therm_TP3$Chl.prop + 1)
W = 0.9211, p-value = 0.003251
##Still not normal but less skewed
##Model with log +1 transformation
#Function of Site and Origin, with Genotype as a Random effect
#Interaction between Site and Origin
Tol_Chl.lme_TP3<-lmer(log(Chl.prop+1)~Origin*Site+(1|Genotype), data=Therm_TP3)
boundary (singular) fit: see help('isSingular')
##Check residuals
Tol_Chl.lme_res_TP3 <- simulateResiduals(fittedModel = Tol_Chl.lme_TP3, plot = F)
plot(Tol_Chl.lme_res_TP3)
##Model results
summary(Tol_Chl.lme_TP3)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(Chl.prop + 1) ~ Origin * Site + (1 | Genotype)
Data: Therm_TP3
REML criterion at convergence: -47.3
Scaled residuals:
Min 1Q Median 3Q Max
-1.7095 -0.6080 -0.1405 0.4999 1.7650
Random effects:
Groups Name Variance Std.Dev.
Genotype (Intercept) 0.00000 0.0000
Residual 0.01594 0.1262
Number of obs: 48, groups: Genotype, 3
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.43851 0.03644 44.00000 12.032 1.65e-15 ***
OriginTransplant 0.03181 0.05154 44.00000 0.617 0.540
SiteSS 0.03852 0.05154 44.00000 0.747 0.459
OriginTransplant:SiteSS -0.02511 0.07289 44.00000 -0.344 0.732
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) OrgnTr SiteSS
OrgnTrnspln -0.707
SiteSS -0.707 0.500
OrgnTrn:SSS 0.500 -0.707 -0.707
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
eta_squared(Tol_Chl.lme_TP3)
# Effect Size for ANOVA (Type III)
Parameter | Eta2 (partial) | 95% CI
-------------------------------------------
Origin | 6.31e-03 | [0.00, 1.00]
Site | 0.01 | [0.00, 1.00]
Origin:Site | 2.69e-03 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_Chl_TP3.res<-data.frame(summary(Tol_Chl.lme_TP3)$coefficients[-1,])
Tol_Chl_TP3.res$Predictor<-c("Origin", "Site", "Origin x Site")
Tol_Chl_TP3.res$EtaSq<-c(eta_squared(Tol_Chl.lme_TP3)$Eta2)
Tol_Chl_TP3.res$Response<-rep("Chlorophyll Retention", nrow(Tol_Chl_TP3.res))
Not normal, but model residuals are OK.
##Pairwise Comparisons of Interest
#Native vs Transplant within a Site
Tol_Chl.pair_TP3<-emmeans(Tol_Chl.lme_TP3, pairwise~Origin | Site)
Tol_Chl.pair_TP3
$emmeans
Site = KL:
Origin emmean SE df lower.CL upper.CL
Native 0.439 0.0364 22.4 0.363 0.514
Transplant 0.470 0.0364 22.4 0.395 0.546
Site = SS:
Origin emmean SE df lower.CL upper.CL
Native 0.477 0.0364 22.4 0.402 0.553
Transplant 0.484 0.0364 22.4 0.408 0.559
Degrees-of-freedom method: kenward-roger
Results are given on the log(mu + 1) (not the response) scale.
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
Native - Transplant -0.03181 0.0515 42 -0.617 0.5404
Site = SS:
contrast estimate SE df t.ratio p.value
Native - Transplant -0.00671 0.0515 42 -0.130 0.8971
Note: contrasts are still on the log(mu + 1) scale
Degrees-of-freedom method: kenward-roger
#Calculate standardized effect sizes for pairwise comparisons
Tol_Chl.pair.es_TP3<-data.frame(eff_size(Tol_Chl.pair_TP3, sigma=sigma(Tol_Chl.lme_TP3), edf=df.residual(Tol_Chl.lme_TP3)))
Since 'object' is a list, we are using the contrasts already present.
Tol_Chl.pair.es_TP3
Tol_Chl.pair.es_TP3<-Tol_Chl.pair.es_TP3 %>% dplyr::rename(contrast.se = contrast, SE.es = SE, df.es = df)
##Save Results
Tol_Chl.pair_TP3.res<-merge(data.frame(Tol_Chl.pair_TP3$contrasts), Tol_Chl.pair.es_TP3[,-c(1)])
Tol_Chl.pair_TP3.res$Response<-rep("Chlorophyll Retention", nrow(Tol_Chl.pair_TP3.res))
##Add Significance levels
Tol_Chl.pair_TP3.res$Sig<-ifelse(Tol_Chl.pair_TP3.res$p.value<0.001, "***",
ifelse(Tol_Chl.pair_TP3.res$p.value<0.01, "**",
ifelse(Tol_Chl.pair_TP3.res$p.value<0.05, "*",
ifelse(Tol_Chl.pair_TP3.res$p.value<0.1, "-", NA))))
Tol_Chl.pair_TP3.res
##Summary statistics by Site and Origin
TP3_TolChl.sum<-summarySE(Therm_TP3, measurevar="Chl.prop", groupvars=c("Site", "Origin", "Site.Orig"), na.rm=TRUE)
##Plot Average Retention across Treatments
TP3_TolChl.plot<-ggplot(TP3_TolChl.sum, aes(x=Site, y=Chl.prop, colour=Site.Orig)) +
scale_colour_manual(values=Orig.colors.o)+
geom_errorbar(aes(ymin=Chl.prop-se, ymax=Chl.prop+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
geom_point(size=point.sz, position=position_dodge(width=0.5))+
theme_classic()+
theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"))+
labs(x="Site and Origin", y="Chlorophyll Retention", colour="Origin")+
ylim(0, 1); TP3_TolChl.plot
#Check normality
hist(Therm_TP3$PAM.prop)
shapiro.test(Therm_TP3$PAM.prop)
Shapiro-Wilk normality test
data: Therm_TP3$PAM.prop
W = 0.96358, p-value = 0.1408
#Normal
##Model
#Function of Site and Origin, with Genotype as a Random effect
#Interaction between Site and Origin
Tol_PAM.lme_TP3<-lmer(PAM.prop~Origin*Site+(1|Genotype), data=Therm_TP3)
##Check residuals
Tol_PAM.lme_res_TP3 <- simulateResiduals(fittedModel = Tol_PAM.lme_TP3, plot = F)
plot(Tol_PAM.lme_res_TP3)
##Model results
summary(Tol_PAM.lme_TP3)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: PAM.prop ~ Origin * Site + (1 | Genotype)
Data: Therm_TP3
REML criterion at convergence: -172.6
Scaled residuals:
Min 1Q Median 3Q Max
-2.29399 -0.71463 0.01706 0.83761 1.77082
Random effects:
Groups Name Variance Std.Dev.
Genotype (Intercept) 3.177e-05 0.005637
Residual 9.061e-04 0.030101
Number of obs: 48, groups: Genotype, 3
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.954325 0.009279 14.520200 102.850 <2e-16 ***
OriginTransplant 0.009025 0.012289 42.000000 0.734 0.467
SiteSS 0.003975 0.012289 42.000000 0.323 0.748
OriginTransplant:SiteSS -0.019133 0.017379 42.000000 -1.101 0.277
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) OrgnTr SiteSS
OrgnTrnspln -0.662
SiteSS -0.662 0.500
OrgnTrn:SSS 0.468 -0.707 -0.707
eta_squared(Tol_PAM.lme_TP3)
# Effect Size for ANOVA (Type III)
Parameter | Eta2 (partial) | 95% CI
-------------------------------------------
Origin | 9.25e-05 | [0.00, 1.00]
Site | 9.76e-03 | [0.00, 1.00]
Origin:Site | 0.03 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_PAM_TP3.res<-data.frame(summary(Tol_PAM.lme_TP3)$coefficients[-1,])
Tol_PAM_TP3.res$Predictor<-c("Origin", "Site", "Origin x Site")
Tol_PAM_TP3.res$EtaSq<-c(eta_squared(Tol_PAM.lme_TP3)$Eta2)
Tol_PAM_TP3.res$Response<-rep("FvFm Retention", nrow(Tol_PAM_TP3.res))
##Pairwise Comparisons of Interest
#Native vs Transplant within a Site
Tol_PAM.pair_TP3<-emmeans(Tol_PAM.lme_TP3, pairwise~Origin | Site)
Tol_PAM.pair_TP3
$emmeans
Site = KL:
Origin emmean SE df lower.CL upper.CL
Native 0.954 0.00928 14.5 0.934 0.974
Transplant 0.963 0.00928 14.5 0.944 0.983
Site = SS:
Origin emmean SE df lower.CL upper.CL
Native 0.958 0.00928 14.5 0.938 0.978
Transplant 0.948 0.00928 14.5 0.928 0.968
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
Native - Transplant -0.00902 0.0123 42 -0.734 0.4668
Site = SS:
contrast estimate SE df t.ratio p.value
Native - Transplant 0.01011 0.0123 42 0.823 0.4154
Degrees-of-freedom method: kenward-roger
#Calculate standardized effect sizes for pairwise comparisons
Tol_PAM.pair.es_TP3<-data.frame(eff_size(Tol_PAM.pair_TP3, sigma=sigma(Tol_PAM.lme_TP3), edf=df.residual(Tol_PAM.lme_TP3)))
Since 'object' is a list, we are using the contrasts already present.
Tol_PAM.pair.es_TP3
Tol_PAM.pair.es_TP3<-Tol_PAM.pair.es_TP3 %>% dplyr::rename(contrast.se = contrast, SE.es = SE, df.es = df)
##Save Results
Tol_PAM.pair_TP3.res<-merge(data.frame(Tol_PAM.pair_TP3$contrasts), Tol_PAM.pair.es_TP3[,-c(1)])
Tol_PAM.pair_TP3.res$Response<-rep("FvFm Retention", nrow(Tol_PAM.pair_TP3.res))
##Add Significance levels
Tol_PAM.pair_TP3.res$Sig<-ifelse(Tol_PAM.pair_TP3.res$p.value<0.001, "***",
ifelse(Tol_PAM.pair_TP3.res$p.value<0.01, "**",
ifelse(Tol_PAM.pair_TP3.res$p.value<0.05, "*",
ifelse(Tol_PAM.pair_TP3.res$p.value<0.1, "-", NA))))
Tol_PAM.pair_TP3.res
##Summary statistics by Site and Origin
TP3_TolPAM.sum<-summarySE(Therm_TP3, measurevar="PAM.prop", groupvars=c("Site", "Origin", "Site.Orig"), na.rm=TRUE)
##Plot Average Retention across Treatments
TP3_TolPAM.plot<-ggplot(TP3_TolPAM.sum, aes(x=Site, y=PAM.prop, colour=Site.Orig)) +
scale_colour_manual(values=Orig.colors.o)+
geom_errorbar(aes(ymin=PAM.prop-se, ymax=PAM.prop+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
geom_point(size=point.sz, position=position_dodge(width=0.5))+
theme_classic()+
theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"))+
labs(x="Site and Origin", y="Photochemical Efficiency Retention", colour="Origin")+
ylim(0.7, 1); TP3_TolPAM.plot
##Combine Results
Tol_TP3.res<-rbind(Tol_PAM_TP3.res, Tol_Sym_TP3.res, Tol_Chl_TP3.res)
Tol.pair.res_TP3<-rbind(Tol_PAM.pair_TP3.res, Tol_Sym.pair_TP3.res, Tol_Chl.pair_TP3.res)
##Add Timepoint
Tol_TP3.res$TimeP<-rep("TP3", nrow(Tol_TP3.res))
Tol.pair.res_TP3$TimeP<-rep("TP3", nrow(Tol.pair.res_TP3))
##Subset Timepoint 4
Therm_TP4<-subset(Therm_H, TimeP=="TP4")
#Check normality
hist(Therm_TP4$Sym.prop)
shapiro.test(Therm_TP4$Sym.prop)
Shapiro-Wilk normality test
data: Therm_TP4$Sym.prop
W = 0.93453, p-value = 0.01005
#Not normal
hist(log(Therm_TP4$Sym.prop+1))
shapiro.test(log(Therm_TP4$Sym.prop+1))
Shapiro-Wilk normality test
data: log(Therm_TP4$Sym.prop + 1)
W = 0.95458, p-value = 0.061
#Normal
##Model with log +1 transformation
#Function of Site and Origin, with Genotype as a Random effect
#Interaction between Site and Origin
Tol_Sym.lme_TP4<-lmer(log(Sym.prop+1)~Origin*Site+(1|Genotype), data=Therm_TP4)
##Check residuals
Tol_Sym.lme_res_TP4 <- simulateResiduals(fittedModel = Tol_Sym.lme_TP4, plot = F)
plot(Tol_Sym.lme_res_TP4)
##Model results
summary(Tol_Sym.lme_TP4)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(Sym.prop + 1) ~ Origin * Site + (1 | Genotype)
Data: Therm_TP4
REML criterion at convergence: -44.7
Scaled residuals:
Min 1Q Median 3Q Max
-1.5412 -0.7855 -0.1534 0.5518 2.2412
Random effects:
Groups Name Variance Std.Dev.
Genotype (Intercept) 0.001295 0.03598
Residual 0.016289 0.12763
Number of obs: 48, groups: Genotype, 3
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.24607 0.04230 9.94436 5.818 0.000173 ***
OriginTransplant 0.04098 0.05210 42.00002 0.787 0.435975
SiteSS -0.03092 0.05210 42.00002 -0.593 0.556032
OriginTransplant:SiteSS -0.03256 0.07369 42.00002 -0.442 0.660816
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) OrgnTr SiteSS
OrgnTrnspln -0.616
SiteSS -0.616 0.500
OrgnTrn:SSS 0.436 -0.707 -0.707
eta_squared(Tol_Sym.lme_TP4)
# Effect Size for ANOVA (Type III)
Parameter | Eta2 (partial) | 95% CI
-------------------------------------------
Origin | 0.01 | [0.00, 1.00]
Site | 0.04 | [0.00, 1.00]
Origin:Site | 4.63e-03 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_Sym_TP4.res<-data.frame(summary(Tol_Sym.lme_TP4)$coefficients[-1,])
Tol_Sym_TP4.res$Predictor<-c("Origin", "Site", "Origin x Site")
Tol_Sym_TP4.res$EtaSq<-c(eta_squared(Tol_Sym.lme_TP4)$Eta2)
Tol_Sym_TP4.res$Response<-rep("Symbiont Retention", nrow(Tol_Sym_TP4.res))
##Pairwise Comparisons of Interest
#Native vs Transplant within a Site
Tol_Sym.pair_TP4<-emmeans(Tol_Sym.lme_TP4, pairwise~Origin | Site)
Tol_Sym.pair_TP4
$emmeans
Site = KL:
Origin emmean SE df lower.CL upper.CL
Native 0.246 0.0423 9.94 0.152 0.340
Transplant 0.287 0.0423 9.94 0.193 0.381
Site = SS:
Origin emmean SE df lower.CL upper.CL
Native 0.215 0.0423 9.94 0.121 0.309
Transplant 0.224 0.0423 9.94 0.129 0.318
Degrees-of-freedom method: kenward-roger
Results are given on the log(mu + 1) (not the response) scale.
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
Native - Transplant -0.04098 0.0521 42 -0.787 0.4360
Site = SS:
contrast estimate SE df t.ratio p.value
Native - Transplant -0.00842 0.0521 42 -0.162 0.8724
Note: contrasts are still on the log(mu + 1) scale
Degrees-of-freedom method: kenward-roger
#Calculate standardized effect sizes for pairwise comparisons
Tol_Sym.pair.es_TP4<-data.frame(eff_size(Tol_Sym.pair_TP4, sigma=sigma(Tol_Sym.lme_TP4), edf=df.residual(Tol_Sym.lme_TP4)))
Since 'object' is a list, we are using the contrasts already present.
Tol_Sym.pair.es_TP4
Tol_Sym.pair.es_TP4<-Tol_Sym.pair.es_TP4 %>% dplyr::rename(contrast.se = contrast, SE.es = SE, df.es = df)
##Save Results
Tol_Sym.pair_TP4.res<-merge(data.frame(Tol_Sym.pair_TP4$contrasts), Tol_Sym.pair.es_TP4[,-c(1)])
Tol_Sym.pair_TP4.res$Response<-rep("Symbiont Retention", nrow(Tol_Sym.pair_TP4.res))
##Add Significance levels
Tol_Sym.pair_TP4.res$Sig<-ifelse(Tol_Sym.pair_TP4.res$p.value<0.001, "***",
ifelse(Tol_Sym.pair_TP4.res$p.value<0.01, "**",
ifelse(Tol_Sym.pair_TP4.res$p.value<0.05, "*",
ifelse(Tol_Sym.pair_TP4.res$p.value<0.1, "-", NA))))
Tol_Sym.pair_TP4.res
##Summary statistics by Site and Origin
TP4_TolSym.sum<-summarySE(Therm_TP4, measurevar="Sym.prop", groupvars=c("Site", "Origin", "Site.Orig"), na.rm=TRUE)
##Plot Average Retention across Treatments
TP4_TolSym.plot<-ggplot(TP4_TolSym.sum, aes(x=Site, y=Sym.prop, colour=Site.Orig)) +
scale_colour_manual(values=Orig.colors.o)+
geom_errorbar(aes(ymin=Sym.prop-se, ymax=Sym.prop+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
geom_point(size=point.sz, position=position_dodge(width=0.5))+
theme_classic()+
theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"))+
labs(x="Site and Origin", y="Symbiont Retention", colour="Origin")+
ylim(0, 1); TP4_TolSym.plot
#Check normality
hist(Therm_TP4$Chl.prop)
shapiro.test(Therm_TP4$Chl.prop)
Shapiro-Wilk normality test
data: Therm_TP4$Chl.prop
W = 0.94204, p-value = 0.02131
#Not normal
hist(log(Therm_TP4$Chl.prop+1))
shapiro.test(log(Therm_TP4$Chl.prop+1))
Shapiro-Wilk normality test
data: log(Therm_TP4$Chl.prop + 1)
W = 0.9512, p-value = 0.04824
##Nearly normal
##Model with log +1 transformation
#Function of Site and Origin, with Genotype as a Random effect
#Interaction between Site and Origin
Tol_Chl.lme_TP4<-lmer(log(Chl.prop+1)~Origin*Site+(1|Genotype), data=Therm_TP4)
boundary (singular) fit: see help('isSingular')
##Check residuals
Tol_Chl.lme_res_TP4 <- simulateResiduals(fittedModel = Tol_Chl.lme_TP4, plot = F)
plot(Tol_Chl.lme_res_TP4)
##Model results
summary(Tol_Chl.lme_TP4)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(Chl.prop + 1) ~ Origin * Site + (1 | Genotype)
Data: Therm_TP4
REML criterion at convergence: -139.8
Scaled residuals:
Min 1Q Median 3Q Max
-1.6397 -0.7647 -0.1492 0.5006 2.4251
Random effects:
Groups Name Variance Std.Dev.
Genotype (Intercept) 0.000000 0.00000
Residual 0.001803 0.04246
Number of obs: 47, groups: Genotype, 3
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.09448 0.01226 43.00000 7.707 1.24e-09 ***
OriginTransplant -0.01266 0.01734 43.00000 -0.730 0.469
SiteSS -0.02659 0.01773 43.00000 -1.500 0.141
OriginTransplant:SiteSS 0.02939 0.02479 43.00000 1.185 0.242
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) OrgnTr SiteSS
OrgnTrnspln -0.707
SiteSS -0.692 0.489
OrgnTrn:SSS 0.494 -0.699 -0.715
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
eta_squared(Tol_Chl.lme_TP4)
# Effect Size for ANOVA (Type III)
Parameter | Eta2 (partial) | 95% CI
-------------------------------------------
Origin | 6.26e-04 | [0.00, 1.00]
Site | 0.02 | [0.00, 1.00]
Origin:Site | 0.03 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_Chl_TP4.res<-data.frame(summary(Tol_Chl.lme_TP4)$coefficients[-1,])
Tol_Chl_TP4.res$Predictor<-c("Origin", "Site", "Origin x Site")
Tol_Chl_TP4.res$EtaSq<-c(eta_squared(Tol_Chl.lme_TP4)$Eta2)
Tol_Chl_TP4.res$Response<-rep("Chlorophyll Retention", nrow(Tol_Chl_TP4.res))
##Pairwise Comparisons of Interest
#Native vs Transplant within a Site
Tol_Chl.pair_TP4<-emmeans(Tol_Chl.lme_TP4, pairwise~Origin | Site)
Tol_Chl.pair_TP4
$emmeans
Site = KL:
Origin emmean SE df lower.CL upper.CL
Native 0.0945 0.0123 21.6 0.0690 0.1199
Transplant 0.0818 0.0123 21.6 0.0564 0.1073
Site = SS:
Origin emmean SE df lower.CL upper.CL
Native 0.0679 0.0129 23.5 0.0413 0.0944
Transplant 0.0846 0.0123 21.6 0.0592 0.1101
Degrees-of-freedom method: kenward-roger
Results are given on the log(mu + 1) (not the response) scale.
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
Native - Transplant 0.0127 0.0173 41.0 0.730 0.4693
Site = SS:
contrast estimate SE df t.ratio p.value
Native - Transplant -0.0167 0.0178 41.2 -0.942 0.3517
Note: contrasts are still on the log(mu + 1) scale
Degrees-of-freedom method: kenward-roger
#Calculate standardized effect sizes for pairwise comparisons
Tol_Chl.pair.es_TP4<-data.frame(eff_size(Tol_Chl.pair_TP4, sigma=sigma(Tol_Chl.lme_TP4), edf=df.residual(Tol_Chl.lme_TP4)))
Since 'object' is a list, we are using the contrasts already present.
Tol_Chl.pair.es_TP4
Tol_Chl.pair.es_TP4<-Tol_Chl.pair.es_TP4 %>% dplyr::rename(contrast.se = contrast, SE.es = SE, df.es = df)
##Save Results
Tol_Chl.pair_TP4.res<-merge(data.frame(Tol_Chl.pair_TP4$contrasts), Tol_Chl.pair.es_TP4[,-c(1)])
Tol_Chl.pair_TP4.res$Response<-rep("Chlorophyll Retention", nrow(Tol_Chl.pair_TP4.res))
##Add Significance levels
Tol_Chl.pair_TP4.res$Sig<-ifelse(Tol_Chl.pair_TP4.res$p.value<0.001, "***",
ifelse(Tol_Chl.pair_TP4.res$p.value<0.01, "**",
ifelse(Tol_Chl.pair_TP4.res$p.value<0.05, "*",
ifelse(Tol_Chl.pair_TP4.res$p.value<0.1, "-", NA))))
Tol_Chl.pair_TP4.res
##Summary statistics by Site and Origin
TP4_TolChl.sum<-summarySE(Therm_TP4, measurevar="Chl.prop", groupvars=c("Site", "Origin", "Site.Orig"), na.rm=TRUE)
##Plot Average Retention across Treatments
TP4_TolChl.plot<-ggplot(TP4_TolChl.sum, aes(x=Site, y=Chl.prop, colour=Site.Orig)) +
scale_colour_manual(values=Orig.colors.o)+
geom_errorbar(aes(ymin=Chl.prop-se, ymax=Chl.prop+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
geom_point(size=point.sz, position=position_dodge(width=0.5))+
theme_classic()+
theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"))+
labs(x="Site and Origin", y="Chlorophyll Retention", colour="Origin")+
ylim(0, 0.35); TP4_TolChl.plot
#Check normality
hist(Therm_TP4$PAM.prop)
shapiro.test(Therm_TP4$PAM.prop)
Shapiro-Wilk normality test
data: Therm_TP4$PAM.prop
W = 0.86379, p-value = 5.159e-05
#Not normal
hist(Therm_TP4$PAM.prop^2)
shapiro.test(Therm_TP4$PAM.prop^2)
Shapiro-Wilk normality test
data: Therm_TP4$PAM.prop^2
W = 0.90996, p-value = 0.00134
##Still not normal
##Model
#Function of Site and Origin, with Genotype as a Random effect
#Interaction between Site and Origin
Tol_PAM.lme_TP4<-lmer(PAM.prop^2~Origin*Site+(1|Genotype), data=Therm_TP4)
##Check residuals
Tol_PAM.lme_res_TP4 <- simulateResiduals(fittedModel = Tol_PAM.lme_TP4, plot = F)
plot(Tol_PAM.lme_res_TP4)
##Model results
summary(Tol_PAM.lme_TP4)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: PAM.prop^2 ~ Origin * Site + (1 | Genotype)
Data: Therm_TP4
REML criterion at convergence: -31
Scaled residuals:
Min 1Q Median 3Q Max
-1.81938 -0.46352 -0.06318 0.64045 2.58461
Random effects:
Groups Name Variance Std.Dev.
Genotype (Intercept) 0.01982 0.1408
Residual 0.02031 0.1425
Number of obs: 48, groups: Genotype, 3
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.58067 0.09109 2.78330 6.374 0.00976 **
OriginTransplant 0.03985 0.05818 42.00000 0.685 0.49718
SiteSS 0.03545 0.05818 42.00000 0.609 0.54566
OriginTransplant:SiteSS -0.07357 0.08228 42.00000 -0.894 0.37633
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) OrgnTr SiteSS
OrgnTrnspln -0.319
SiteSS -0.319 0.500
OrgnTrn:SSS 0.226 -0.707 -0.707
eta_squared(Tol_PAM.lme_TP4)
# Effect Size for ANOVA (Type III)
Parameter | Eta2 (partial) | 95% CI
-------------------------------------------
Origin | 1.32e-04 | [0.00, 1.00]
Site | 2.53e-05 | [0.00, 1.00]
Origin:Site | 0.02 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_PAM_TP4.res<-data.frame(summary(Tol_PAM.lme_TP4)$coefficients[-1,])
Tol_PAM_TP4.res$Predictor<-c("Origin", "Site", "Origin x Site")
Tol_PAM_TP4.res$EtaSq<-c(eta_squared(Tol_PAM.lme_TP4)$Eta2)
Tol_PAM_TP4.res$Response<-rep("FvFm Retention", nrow(Tol_PAM_TP4.res))
##Pairwise Comparisons of Interest
#Native vs Transplant within a Site
Tol_PAM.pair_TP4<-emmeans(Tol_PAM.lme_TP4, pairwise~Origin | Site)
Tol_PAM.pair_TP4
$emmeans
Site = KL:
Origin emmean SE df lower.CL upper.CL
Native 0.581 0.0911 2.78 0.278 0.884
Transplant 0.621 0.0911 2.78 0.317 0.924
Site = SS:
Origin emmean SE df lower.CL upper.CL
Native 0.616 0.0911 2.78 0.313 0.919
Transplant 0.582 0.0911 2.78 0.279 0.885
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
Native - Transplant -0.0398 0.0582 42 -0.685 0.4972
Site = SS:
contrast estimate SE df t.ratio p.value
Native - Transplant 0.0337 0.0582 42 0.580 0.5652
Degrees-of-freedom method: kenward-roger
#Calculate standardized effect sizes for pairwise comparisons
Tol_PAM.pair.es_TP4<-data.frame(eff_size(Tol_PAM.pair_TP4, sigma=sigma(Tol_PAM.lme_TP4), edf=df.residual(Tol_PAM.lme_TP4)))
Since 'object' is a list, we are using the contrasts already present.
Tol_PAM.pair.es_TP4
Tol_PAM.pair.es_TP4<-Tol_PAM.pair.es_TP4 %>% dplyr::rename(contrast.se = contrast, SE.es = SE, df.es = df)
##Save Results
Tol_PAM.pair_TP4.res<-merge(data.frame(Tol_PAM.pair_TP4$contrasts), Tol_PAM.pair.es_TP4[,-c(1)])
Tol_PAM.pair_TP4.res$Response<-rep("FvFm Retention", nrow(Tol_PAM.pair_TP4.res))
##Add Significance levels
Tol_PAM.pair_TP4.res$Sig<-ifelse(Tol_PAM.pair_TP4.res$p.value<0.001, "***",
ifelse(Tol_PAM.pair_TP4.res$p.value<0.01, "**",
ifelse(Tol_PAM.pair_TP4.res$p.value<0.05, "*",
ifelse(Tol_PAM.pair_TP4.res$p.value<0.1, "-", NA))))
Tol_PAM.pair_TP4.res
##Summary statistics by Site and Origin
TP4_TolPAM.sum<-summarySE(Therm_TP4, measurevar="PAM.prop", groupvars=c("Site", "Origin", "Site.Orig"), na.rm=TRUE)
##Plot Average Retention across Treatments
TP4_TolPAM.plot<-ggplot(TP4_TolPAM.sum, aes(x=Site, y=PAM.prop, colour=Site.Orig)) +
scale_colour_manual(values=Orig.colors.o)+
geom_errorbar(aes(ymin=PAM.prop-se, ymax=PAM.prop+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
geom_point(size=point.sz, position=position_dodge(width=0.5))+
theme_classic()+
theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"))+
labs(x="Site and Origin", y="Photochemical Efficiency Retention", colour="Origin")+
ylim(0.6, 1); TP4_TolPAM.plot
##Combine Results
Tol_TP4.res<-rbind(Tol_PAM_TP4.res, Tol_Sym_TP4.res, Tol_Chl_TP4.res)
Tol.pair.res_TP4<-rbind(Tol_PAM.pair_TP4.res, Tol_Sym.pair_TP4.res, Tol_Chl.pair_TP4.res)
##Add Timepoint
Tol_TP4.res$TimeP<-rep("TP4", nrow(Tol_TP4.res))
Tol.pair.res_TP4$TimeP<-rep("TP4", nrow(Tol.pair.res_TP4))
Focus on Fv/Fm at KL TP 1 and 2 where difference were significant
##Summary statistics by Genotype, Origin, Site, and Timepoint
Tol.PAM.sum<-summarySE(Therm_H, measurevar="PAM.prop", groupvars=c("Genotype", "Origin", "Site", "TimeP"), na.rm=TRUE)
##Remove Initial Timepoint
Tol.PAM.sum<-Tol.PAM.sum[-c(which(Tol.PAM.sum$TimeP=="IN")),]
##Calculate Percent Difference between Native and Transplant by Genotype
Tol.PAM.dif<-Tol.PAM.sum[which(Tol.PAM.sum$Origin=="Native"), c(1,3,4,6)]
Tol.PAM.dif<-Tol.PAM.dif %>% dplyr::rename(PAM_N = PAM.prop)
Tol.PAM.dif$PAM_T<-Tol.PAM.sum$PAM.prop[which(Tol.PAM.sum$Origin=="Transplant")]
##KL
Tol.PAM.dif_KL<-subset(Tol.PAM.dif, Site=="KL")
Tol.PAM.dif_KL$Perc.Inc<-((Tol.PAM.dif_KL$PAM_T-Tol.PAM.dif_KL$PAM_N)/Tol.PAM.dif_KL$PAM_N)*100
Tol.PAM.dif_KL
##Mean and Standard Error by Timepoint
Tol.PAM.dif_KL.sum<-summarySE(Tol.PAM.dif_KL, measurevar="Perc.Inc", groupvars=c("TimeP"), na.rm=TRUE)
Tol.PAM.dif_KL.sum
##Combine Results
Growth.TP.res<-rbind(TLE_TP1.2.res, TLE_TP3.4.res)
Growth.TP.res<-Growth.TP.res %>% dplyr::rename(TimeP = Timepoint)
Tol_TP.res<-rbind(Tol_TP1.res, Tol_TP2.res, Tol_TP3.res, Tol_TP4.res)
Perf.res<-rbind(Growth.TP.res, Tol_TP.res)
Perf.res<-Perf.res %>% dplyr::rename(p = Pr...t..)
##Responses for Figure
Perf.res.plot<-Perf.res[-c(which(Perf.res$Response=="Chlorophyll Retention")),]
##Add Significance levels
Perf.res.plot$Sig<-ifelse(Perf.res.plot$p<0.001, "***", ifelse(Perf.res.plot$p<0.01, "**",
ifelse(Perf.res.plot$p<0.05, "*", ifelse(Perf.res.plot$p<0.1, "-", NA))))
##Set Order
Perf.res.plot$Response<-factor(Perf.res.plot$Response, levels=c("Growth", "FvFm Retention", "Symbiont Retention"), ordered=TRUE)
##Match Growth Timepoints
Perf.res.plot$TimeP[which(Perf.res.plot$TimeP=="TP1_2")]<-"TP2"
Perf.res.plot$TimeP[which(Perf.res.plot$TimeP=="TP3_4")]<-"TP4"
Perf.Site.Orig.ES.plot<-ggplot(Perf.res.plot, aes(fill=Predictor, y=EtaSq, x=TimeP)) +
geom_bar(position="stack", stat="identity")+
facet_wrap(vars(Response))+
theme_bw()+
scale_fill_manual(values=c("indianred2", "mediumpurple3", "steelblue2"))+
theme(strip.text=element_text(size=axis.title.sz),
strip.background = element_rect(fill="white"),
axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"),
legend.position=c(.08, .8),
panel.grid.major.x=element_blank())+
labs(x="Time Point", y=expression(paste("Predictor Effect Size (p", eta^2, ")")), colour="Predictor")+
geom_text(aes(label=Sig), position=position_stack(vjust=0.5), size=sig.sz); Perf.Site.Orig.ES.plot
##Combine Pairwise Results
Growth.pair.res<-rbind(TLE.pair_TP1.2.res, TLE.pair_TP3.4.res)
Tol.pair.res<-rbind(Tol.pair.res_TP1, Tol.pair.res_TP2, Tol.pair.res_TP3, Tol.pair.res_TP4)
Perf.pair.res<-rbind(Growth.pair.res[,c(1:13, 15, 14)], Tol.pair.res)
##Responses for Figure
Perf.pair.res.plot<-Perf.pair.res[-c(which(Perf.pair.res$Response=="Chlorophyll Retention")),]
##Add Significance
Perf.pair.res.plot$Significant<-ifelse(Perf.pair.res.plot$p.value<0.05, "sig", "non")
##Set Order
Perf.pair.res.plot$Response<-factor(Perf.pair.res.plot$Response, levels=c("Growth", "FvFm Retention", "Symbiont Retention"), ordered=TRUE)
##Match Growth Timepoints
Perf.pair.res.plot$TimeP[which(Perf.pair.res.plot$TimeP=="TP1_2")]<-"TP2"
Perf.pair.res.plot$TimeP[which(Perf.pair.res.plot$TimeP=="TP3_4")]<-"TP4"
##Set Factors
Perf.pair.res.plot$Site<-factor(Perf.pair.res.plot$Site, levels=c("KL", "SS"))
##Plot Effect Size across Timepoints
Perf.Orig.ES.plot<-ggplot(Perf.pair.res.plot, aes(x=TimeP, y=effect.size)) +
geom_rect(fill="grey90", alpha=0.4, xmin=-Inf, xmax=Inf, ymin=-.4, ymax=0.4)+
geom_rect(fill="grey80", alpha=0.4, xmin=-Inf, xmax=Inf, ymin=-.25, ymax=0.25)+
geom_rect(fill="grey70", alpha=0.4, xmin=-Inf, xmax=Inf, ymin=-.1, ymax=0.1)+
geom_hline(yintercept=0, linetype="dashed", color="grey50", linewidth=0.5)+
geom_errorbar(aes(ymin=effect.size-SE.es, ymax=effect.size+SE.es, colour=Site), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
geom_point(aes(shape=Significant, colour=Site), size=point.sz, position=position_dodge(width=0.5))+
facet_wrap(vars(Response))+
scale_colour_manual(values=Site.colors.o)+
scale_shape_manual(values=c(1, 16))+
theme_bw()+
theme(axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"),
legend.position=c(.1, .1),
legend.direction = "horizontal",
strip.text=element_text(size=axis.title.sz),
strip.background = element_rect(fill="white"),
panel.grid.major.y=element_blank())+
labs(x="Time Point", y="Origin Effect Size (Cohen's d)", colour="Site")+
guides(shape="none"); Perf.Orig.ES.plot
Need to fix dodge for points and error bars
Calculating the average of Growth (TP1 to TP2) and Thermal Tolerance (Retention of Symbionts, Chlorophyll, and Photochemical Efficiency at TP1) for each Set (Site, Genotype, Origin).
##Average
names(Growth_TP1.2)
[1] "Site" "ID" "Orig" "Origin" "Genotype" "TL.TP1_cm" "TL.TP2_cm"
[8] "TL.TP3_cm" "TL.TP4_cm" "Set" "Site.Orig" "TP1" "TP2" "TP1.2_days"
[15] "Ext_cm" "TLE_cm.day"
TP1.2_Growth_Set<-aggregate(Growth_TP1.2$TLE_cm.day, list(Growth_TP1.2$Site, Growth_TP1.2$Genotype, Growth_TP1.2$Orig), mean)
names(TP1.2_Growth_Set)<-c("Site", "Genotype", "Orig", "TLE_cm.day")
##Add Standard Error
TP1.2_Growth_Set$TLE_se<-aggregate(Growth_TP1.2$TLE_cm.day, list(Growth_TP1.2$Site, Growth_TP1.2$Genotype, Growth_TP1.2$Orig), std.error)[,4]
##Average
names(Therm_TP1)
[1] "TimeP" "Site" "Genotype" "Origin" "ID"
[6] "RandN" "Treat" "Treatment" "Orig" "Set"
[11] "Site.Orig" "SA_cm2" "Chl_ug.cm2" "Sym10.6_cm2" "Fv_Fm"
[16] "Chl_ug.cm2_C" "Sym10.6_cm2_C" "Fv_Fm_C" "Chl.prop" "Sym.prop"
[21] "PAM.prop"
T1_Therm_Set<-aggregate(Therm_TP1[,c(19:21)], list(Therm_TP1$Site, Therm_TP1$Genotype, Therm_TP1$Orig), mean)
names(T1_Therm_Set)<-c("Site", "Genotype", "Orig", "Chl.prop", "Sym.prop", "PAM.prop")
##Add Standard Error
T1_Therm_Set$Chl_se<-aggregate(Therm_TP1$Chl.prop, list(Therm_TP1$Site, Therm_TP1$Genotype, Therm_TP1$Orig), std.error)[,4]
T1_Therm_Set$Sym_se<-aggregate(Therm_TP1$Sym.prop, list(Therm_TP1$Site, Therm_TP1$Genotype, Therm_TP1$Orig), std.error)[,4]
T1_Therm_Set$PAM_se<-aggregate(Therm_TP1$PAM.prop, list(Therm_TP1$Site, Therm_TP1$Genotype, Therm_TP1$Orig), std.error)[,4]
TP1.2_Perf_Set<-merge(TP1.2_Growth_Set, T1_Therm_Set, all=TRUE)
TP1.2_Perf_Set$Set<-paste(TP1.2_Perf_Set$Site, TP1.2_Perf_Set$Genotype, TP1.2_Perf_Set$Orig, sep=".")
##Set factor variables
TP1.2_Perf_Set$Site<-factor(TP1.2_Perf_Set$Site, levels=c("KL", "SS"))
TP1.2_Perf_Set$Genotype<-factor(TP1.2_Perf_Set$Genotype, levels=c("AC8", "AC10", "AC12"))
TP1.2_Perf_Set$Orig<-factor(TP1.2_Perf_Set$Orig, levels=c("N", "T"))
##Add a Sample Set Variable
TP1.2_Perf_Set$Set<-paste(TP1.2_Perf_Set$Site, TP1.2_Perf_Set$Genotype, TP1.2_Perf_Set$Orig, sep=".")
TP1.2_Perf_Set$Set<-factor(TP1.2_Perf_Set$Set)
##Add Site.Orig variable
TP1.2_Perf_Set$Site.Orig<-paste(TP1.2_Perf_Set$Site, TP1.2_Perf_Set$Orig, sep=".")
TP1.2_Perf_Set$Site.Orig<-factor(TP1.2_Perf_Set$Site.Orig, levels=c("KL.N", "KL.T", "SS.N", "SS.T"))
cor.test(TP1.2_Perf_Set$TLE_cm.day, TP1.2_Perf_Set$Sym.prop, method="spearman")
Spearman's rank correlation rho
data: TP1.2_Perf_Set$TLE_cm.day and TP1.2_Perf_Set$Sym.prop
S = 462, p-value = 0.03733
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
-0.6153846
TP1.2_TLE_Sym.plot<-ggplot(TP1.2_Perf_Set, aes(x=TLE_cm.day, y=Sym.prop)) +
geom_smooth(method="lm", color="#AA185AFF", fill="#F8D7BFFF", se=TRUE)+
geom_point(aes(colour=Site.Orig, shape=Genotype),size=point.sz)+
scale_colour_manual(values=Orig.colors.o)+
scale_shape_manual(values=Geno.shapes.o)+
theme_classic()+
theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"))+
labs(x=expression(paste('Growth Rate (cm day'^-1*")")),
y="Symbiont Retention",
colour="Site Origin")+
annotate("text", x = 0.07, y = -.15, label=expression(bolditalic(paste(r[S], " = -0.615, p = 0.037"))), size=sig.sz, hjust = 0); TP1.2_TLE_Sym.plot
cor.test(TP1.2_Perf_Set$TLE_cm.day, TP1.2_Perf_Set$Chl.prop, method="spearman")
Spearman's rank correlation rho
data: TP1.2_Perf_Set$TLE_cm.day and TP1.2_Perf_Set$Chl.prop
S = 454, p-value = 0.04884
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
-0.5874126
TP1.2_TLE_Chl.plot<-ggplot(TP1.2_Perf_Set, aes(x=TLE_cm.day, y=Chl.prop)) +
geom_smooth(method="lm", color="#AA185AFF", fill="#F8D7BFFF", se=TRUE)+
geom_point(aes(colour=Site.Orig, shape=Genotype),size=point.sz)+
scale_colour_manual(values=Orig.colors.o)+
scale_shape_manual(values=Geno.shapes.o)+
theme_classic()+
theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"))+
labs(x=expression(paste('Growth Rate (cm day'^-1*")")),
y="Chlorophyll Retention",
colour="Site Origin")+
annotate("text", x = 0.07, y = -.05, label=expression(bolditalic(paste(r[S], " = -0.587, p = 0.049"))), size=sig.sz, hjust = 0); TP1.2_TLE_Chl.plot
cor.test(TP1.2_Perf_Set$TLE_cm.day, TP1.2_Perf_Set$PAM.prop, method="spearman")
Spearman's rank correlation rho
data: TP1.2_Perf_Set$TLE_cm.day and TP1.2_Perf_Set$PAM.prop
S = 480, p-value = 0.01883
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
-0.6783217
TP1.2_TLE_PAM.plot<-ggplot(TP1.2_Perf_Set, aes(x=TLE_cm.day, y=PAM.prop)) +
geom_smooth(method="lm", color="#AA185AFF", fill="#F8D7BFFF", se=TRUE)+
geom_point(aes(colour=Site.Orig, shape=Genotype),size=point.sz)+
scale_colour_manual(values=Orig.colors.o)+
scale_shape_manual(values=Geno.shapes.o)+
theme_classic()+
theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"))+
labs(x=expression(paste('Growth Rate (cm day'^-1*")")),
y="Photo. Effic. Retention",
colour="Site Origin")+
annotate("text", x = 0.07, y = 0.7, label=expression(bolditalic(paste(r[S], " = -0.678, p = 0.019"))), size=sig.sz, hjust = 0); TP1.2_TLE_PAM.plot
Calculating the average of Growth (TP3 to TP4) and Thermal Tolerance (Retention of Symbionts, Chlorophyll, and Photochemical Efficiency at TP3) for each Set (Site, Genotype, Origin).
##Average
names(Growth_TP3.4)
[1] "Site" "ID" "Orig" "Origin" "Genotype" "TL.TP1_cm" "TL.TP2_cm"
[8] "TL.TP3_cm" "TL.TP4_cm" "Set" "Site.Orig" "TP3" "TP4" "TP3.4_days"
[15] "Ext_cm" "TLE_cm.day"
TP3.4_Growth_Set<-aggregate(Growth_TP3.4$TLE_cm.day, list(Growth_TP3.4$Site, Growth_TP3.4$Genotype, Growth_TP3.4$Orig), mean)
names(TP3.4_Growth_Set)<-c("Site", "Genotype", "Orig", "TLE_cm.day")
##Add Standard Error
TP3.4_Growth_Set$TLE_se<-aggregate(Growth_TP3.4$TLE_cm.day, list(Growth_TP3.4$Site, Growth_TP3.4$Genotype, Growth_TP3.4$Orig), std.error)[,4]
##Average
names(Therm_TP4)
[1] "TimeP" "Site" "Genotype" "Origin" "ID"
[6] "RandN" "Treat" "Treatment" "Orig" "Set"
[11] "Site.Orig" "SA_cm2" "Chl_ug.cm2" "Sym10.6_cm2" "Fv_Fm"
[16] "Chl_ug.cm2_C" "Sym10.6_cm2_C" "Fv_Fm_C" "Chl.prop" "Sym.prop"
[21] "PAM.prop"
T4_Therm_Set<-aggregate(Therm_TP4[,c(19:21)], list(Therm_TP4$Site, Therm_TP4$Genotype, Therm_TP4$Orig), mean)
names(T4_Therm_Set)<-c("Site", "Genotype", "Orig", "Chl.prop", "Sym.prop", "PAM.prop")
##Add Standard Error
T4_Therm_Set$Chl_se<-aggregate(Therm_TP4$Chl.prop, list(Therm_TP4$Site, Therm_TP4$Genotype, Therm_TP4$Orig), std.error)[,4]
T4_Therm_Set$Sym_se<-aggregate(Therm_TP4$Sym.prop, list(Therm_TP4$Site, Therm_TP4$Genotype, Therm_TP4$Orig), std.error)[,4]
T4_Therm_Set$PAM_se<-aggregate(Therm_TP4$PAM.prop, list(Therm_TP4$Site, Therm_TP4$Genotype, Therm_TP4$Orig), std.error)[,4]
TP3.4_Perf_Set<-merge(TP3.4_Growth_Set, T4_Therm_Set, all=TRUE)
TP3.4_Perf_Set$Set<-paste(TP3.4_Perf_Set$Site, TP3.4_Perf_Set$Genotype, TP3.4_Perf_Set$Orig, sep=".")
##Set factor variables
TP3.4_Perf_Set$Site<-factor(TP3.4_Perf_Set$Site, levels=c("KL", "SS"))
TP3.4_Perf_Set$Genotype<-factor(TP3.4_Perf_Set$Genotype, levels=c("AC8", "AC10", "AC12"))
TP3.4_Perf_Set$Orig<-factor(TP3.4_Perf_Set$Orig, levels=c("N", "T"))
##Add a Sample Set Variable
TP3.4_Perf_Set$Set<-paste(TP3.4_Perf_Set$Site, TP3.4_Perf_Set$Genotype, TP3.4_Perf_Set$Orig, sep=".")
TP3.4_Perf_Set$Set<-factor(TP3.4_Perf_Set$Set)
##Add Site.Orig variable
TP3.4_Perf_Set$Site.Orig<-paste(TP3.4_Perf_Set$Site, TP3.4_Perf_Set$Orig, sep=".")
TP3.4_Perf_Set$Site.Orig<-factor(TP3.4_Perf_Set$Site.Orig, levels=c("KL.N", "KL.T", "SS.N", "SS.T"))
cor.test(TP3.4_Perf_Set$TLE_cm.day, TP3.4_Perf_Set$Sym.prop, method="spearman")
Spearman's rank correlation rho
data: TP3.4_Perf_Set$TLE_cm.day and TP3.4_Perf_Set$Sym.prop
S = 444, p-value = 0.06663
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
-0.5524476
TP3.4_TLE_Sym.plot<-ggplot(TP3.4_Perf_Set, aes(x=TLE_cm.day, y=Sym.prop)) +
geom_smooth(method="lm", color="#AA185AFF", fill="#F8D7BFFF", se=TRUE)+
geom_point(aes(colour=Site.Orig, shape=Genotype),size=point.sz)+
scale_colour_manual(values=Orig.colors.o)+
scale_shape_manual(values=Geno.shapes.o)+
theme_classic()+
theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"))+
labs(x=expression(paste('Growth Rate (cm day'^-1*")")),
y="Symbiont Retention",
colour="Site Origin")+
annotate("text", x = 0.07, y = 0.1, label=expression(italic(paste(r[S], " = -0.552, p = 0.067"))), size=sig.sz, hjust = 0); TP3.4_TLE_Sym.plot
cor.test(TP3.4_Perf_Set$TLE_cm.day, TP3.4_Perf_Set$Chl.prop, method="spearman")
Spearman's rank correlation rho
data: TP3.4_Perf_Set$TLE_cm.day and TP3.4_Perf_Set$Chl.prop
S = 148, p-value = 0.327
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
0.3272727
TP3.4_TLE_Chl.plot<-ggplot(TP3.4_Perf_Set, aes(x=TLE_cm.day, y=Chl.prop)) +
geom_smooth(method="lm", color="#AA185AFF", fill="#F8D7BFFF", se=TRUE)+
geom_point(aes(colour=Site.Orig, shape=Genotype),size=point.sz)+
scale_colour_manual(values=Orig.colors.o)+
scale_shape_manual(values=Geno.shapes.o)+
theme_classic()+
theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"))+
labs(x=expression(paste('Growth Rate (cm day'^-1*")")),
y="Chlorophyll Retention",
colour="Site Origin")+
annotate("text", x = 0.07, y = 0.065, label=expression(italic(paste(r[S], " = 0.327, p = 0.327"))), size=sig.sz, hjust = 0); TP3.4_TLE_Chl.plot
cor.test(TP3.4_Perf_Set$TLE_cm.day, TP3.4_Perf_Set$PAM.prop, method="spearman")
Spearman's rank correlation rho
data: TP3.4_Perf_Set$TLE_cm.day and TP3.4_Perf_Set$PAM.prop
S = 500, p-value = 0.007353
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
-0.7482517
TP3.4_TLE_PAM.plot<-ggplot(TP3.4_Perf_Set, aes(x=TLE_cm.day, y=PAM.prop)) +
geom_smooth(method="lm", color="#AA185AFF", fill="#F8D7BFFF", se=TRUE)+
geom_point(aes(colour=Site.Orig, shape=Genotype),size=point.sz)+
scale_colour_manual(values=Orig.colors.o)+
scale_shape_manual(values=Geno.shapes.o)+
theme_classic()+
theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"))+
labs(x=expression(paste('Growth Rate (cm day'^-1*")")),
y="Photo. Effic. Retention",
colour="Site Origin")+
annotate("text", x = 0.07, y = 0.5, label=expression(bolditalic(paste(r[S], " = -0.748, p = 0.007"))), size=sig.sz, hjust = 0); TP3.4_TLE_PAM.plot
Perf.Site.Orig.ES.plot<-Perf.Site.Orig.ES.plot+ labs(x="")
##Create Panel
Perf.ES_fig<-plot_grid(Perf.Site.Orig.ES.plot, Perf.Orig.ES.plot,
rel_widths=1, rel_heights = 1,
nrow=2, ncol=1, byrow=T, labels = c("A", "B"),
label_size=panel.lab.sz)
Warning: Removed 18 rows containing missing values (`geom_text()`).
##Save Figure
ggsave(filename="Figures/03_Performance/Fig4_PerformanceEffectSizes.png", plot=Perf.ES_fig, dpi=300, width=12, height=9, units="in")
##Timepoint 1 to 2
TP1.2_TLE_PAM.plot<-TP1.2_TLE_PAM.plot + ggtitle("Photo. Effic.")+ labs(x=NULL, y = "Tolerance (Retention)")+
theme(plot.title = element_text(colour="black", size=panel.lab.sz, hjust=0.5), legend.position="none")
TP1.2_TLE_Sym.plot<- TP1.2_TLE_Sym.plot + ggtitle("Symbionts")+ labs(x=NULL, y = NULL)+
theme(plot.title = element_text(colour="black", size=panel.lab.sz, hjust=0.5)) + ylim(-.15, .8)
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
##Timepoint 3 to 4
TP3.4_TLE_PAM.plot<-TP3.4_TLE_PAM.plot + labs(y = "Tolerance (Retention)") +
theme(legend.position="none")
TP3.4_TLE_Sym.plot<-TP3.4_TLE_Sym.plot + labs(y = NULL)
##Create Panel
TradeOff_fig<-plot_grid(TP1.2_TLE_PAM.plot, TP1.2_TLE_Sym.plot,
TP3.4_TLE_PAM.plot, TP3.4_TLE_Sym.plot,
rel_widths=c(.75, 1), rel_heights = c(1, .98),
nrow=2, ncol=2, byrow=T,
labels = c("A", "B", "C", "D"),
label_size=panel.lab.sz)
`geom_smooth()` using formula = 'y ~ x'Warning: is.na() applied to non-(list or vector) of type 'expression'`geom_smooth()` using formula = 'y ~ x'Warning: Removed 1 rows containing non-finite values (`stat_smooth()`).Warning: Removed 1 rows containing missing values (`geom_point()`).Warning: is.na() applied to non-(list or vector) of type 'expression'`geom_smooth()` using formula = 'y ~ x'Warning: is.na() applied to non-(list or vector) of type 'expression'`geom_smooth()` using formula = 'y ~ x'Warning: is.na() applied to non-(list or vector) of type 'expression'
##Save Figure
ggsave(filename="Figures/03_Performance/Fig5_TradeOff.png", plot=TradeOff_fig, dpi=300, width=8, height=8, units="in")
##TP 1 to TP 2
TP1.2_Growth.plot<-TP1.2_Growth.plot+ ggtitle("TP1 to TP2")+
theme(plot.title = element_text(colour="black",
size=panel.lab.sz, hjust=0.5),
legend.position=c(0.25, 0.85), legend.direction="horizontal")+
labs(colour=NULL)+
guides(color=guide_legend(nrow=2, byrow=FALSE))
##TP 3 to TP 4
TP3.4_Growth.plot<-TP3.4_Growth.plot+ ggtitle("TP3 to TP4")+
theme(plot.title = element_text(colour="black",
size=panel.lab.sz, hjust=0.5),
legend.position="none")
##Create Panel
Growth_fig<-plot_grid(TP1.2_Growth.plot, TP3.4_Growth.plot,
rel_widths=1, rel_heights = 1,
nrow=1, ncol=2, byrow=T, labels = c("A", "B"),
label_size=panel.lab.sz)
##Save Figure
ggsave(filename="Figures/03_Performance/FigS5_Growth.png", plot=Growth_fig, dpi=300, width=10, height=4.5, units="in")
##FvFm
IN_PAM.plot<-IN_PAM.plot + labs(x="", y="Fv/Fm") + theme(legend.position="none") + ylim(0.5, 0.655)
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
IN_TolPAM.plot<-IN_TolPAM.plot + labs(x="", y="Retention") +
theme(legend.position="none") + ylim(0.75, 1)
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
##Symbionts
IN_Sym.plot<-IN_Sym.plot + labs(x="", y="Symbionts") + theme(legend.position="none") + ylim(0,1)
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
IN_TolSym.plot<-IN_TolSym.plot + labs(x="", y="Retention") +
theme(legend.position="none") + ylim(.25, .75)
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
##Chlorophyll
IN_Chl.plot<-IN_Chl.plot + labs(x="Site and Treatment", y="Chlorophyll") +
theme(legend.position="none") + ylim(0,1.25)
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
IN_TolChl.plot<-IN_TolChl.plot + labs(x="Site", y="Retention") +
theme(legend.position="none") + ylim(0, 0.25)
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
##Create Panel
IN_HeatAssay_fig<-plot_grid(IN_PAM.plot, IN_TolPAM.plot,
IN_Sym.plot, IN_TolSym.plot,
IN_Chl.plot, IN_TolChl.plot,
rel_widths=c(1, 0.75),
rel_heights = c(1, 1, 1),
nrow=3, ncol=2, byrow=T,
labels = c("A", "B", "C", "D", "E", "F"),
label_size=panel.lab.sz)
##Save Figure
ggsave(filename="Figures/03_Performance/FigS6_InitialHeatAssay.png", plot=IN_HeatAssay_fig, dpi=300, width=8, height=11, units="in")
##FvFm Retention
TP1_TolPAM.plot<-TP1_TolPAM.plot + ggtitle("TP1")+
labs(x=NULL, y="Photo. Effic. Retention", colour=NULL)+
theme(plot.title = element_text(colour="black", size=panel.lab.sz, face="bold", hjust=0.5),
legend.position=c(0.5, 0.15), legend.direction="horizontal")+
guides(color=guide_legend(nrow=2, byrow=FALSE))
TP2_TolPAM.plot<-TP2_TolPAM.plot+ ggtitle("TP2")+ labs(x=NULL, y=NULL) +
theme(plot.title = element_text(colour="black", size=panel.lab.sz, face="bold", hjust=0.5),
legend.position="none")
TP3_TolPAM.plot<-TP3_TolPAM.plot+ ggtitle("TP3")+ labs(x=NULL, y=NULL) +
theme(plot.title = element_text(colour="black", size=panel.lab.sz, face="bold", hjust=0.5),
legend.position="none")
TP4_TolPAM.plot<-TP4_TolPAM.plot+ ggtitle("TP4")+ labs(x=NULL, y=NULL) +
theme(plot.title = element_text(colour="black", size=panel.lab.sz, face="bold", hjust=0.5),
legend.position="none")
##Symbiont Retention
TP1_TolSym.plot<- TP1_TolSym.plot + labs(x=NULL) + theme(legend.position="none")
TP2_TolSym.plot<- TP2_TolSym.plot + labs(x=NULL, y=NULL) + theme(legend.position="none")
TP3_TolSym.plot<- TP3_TolSym.plot + labs(x=NULL, y=NULL) + theme(legend.position="none")
TP4_TolSym.plot<- TP4_TolSym.plot + labs(x=NULL, y=NULL) + theme(legend.position="none")
##Chlorophyll Retention
TP1_TolChl.plot<-TP1_TolChl.plot + theme(legend.position="none")
TP2_TolChl.plot<-TP2_TolChl.plot + labs(y=NULL) + theme(legend.position="none")
TP3_TolChl.plot<-TP3_TolChl.plot + labs(y=NULL) + theme(legend.position="none")
TP4_TolChl.plot<-TP4_TolChl.plot + labs(y=NULL) + theme(legend.position="none")
##Create Panel
Tolerance_fig<-plot_grid(TP1_TolPAM.plot, TP2_TolPAM.plot,
TP3_TolPAM.plot, TP4_TolPAM.plot,
TP1_TolSym.plot, TP2_TolSym.plot,
TP3_TolSym.plot, TP4_TolSym.plot,
TP1_TolChl.plot, TP2_TolChl.plot,
TP3_TolChl.plot, TP4_TolChl.plot,
rel_widths=1, rel_heights = 1,
nrow=3, ncol=4, byrow=T, labels = NULL)
##Save Figure
ggsave(filename="Figures/03_Performance/FigS7_Tolerance.png", plot=Tolerance_fig, dpi=300, width=14, height=12, units="in")
##Combine Results Tables
TableS6_Perf.LM<-Perf.res
##Organize
names(TableS6_Perf.LM)
[1] "Estimate" "Std..Error" "df" "t.value" "p" "Predictor" "EtaSq"
[8] "Response" "TimeP"
TableS6_Perf.LM<-TableS6_Perf.LM %>% dplyr::rename(SE = Std..Error, DF = df, t = t.value, Timepoint = TimeP)
TableS6_Perf.LM<-TableS6_Perf.LM[,c("Timepoint", "Response", "Predictor", "Estimate", "SE", "DF", "t", "p", "EtaSq")]
#Round to 3 digits
TableS6_Perf.LM$Estimate<-round(TableS6_Perf.LM$Estimate, 3)
TableS6_Perf.LM$SE<-round(TableS6_Perf.LM$SE, 3)
TableS6_Perf.LM$t<-round(TableS6_Perf.LM$t, 3)
TableS6_Perf.LM$p<-round(TableS6_Perf.LM$p, 3)
TableS6_Perf.LM$EtaSq<-round(TableS6_Perf.LM$EtaSq, 3)
#Integer
TableS6_Perf.LM$DF<-round(TableS6_Perf.LM$DF, 0)
##Write Out Table
write.csv(TableS6_Perf.LM, "Tables/TableS6_Performance_LM_Results.csv", row.names=FALSE)
##Save Results Table
TableS7_Perf.pair<-Perf.pair.res
##Organize
names(TableS7_Perf.pair)
[1] "Site" "contrast" "estimate" "SE" "df" "t.ratio"
[7] "p.value" "effect.size" "SE.es" "df.es" "lower.CL" "upper.CL"
[13] "Response" "Sig" "TimeP"
TableS7_Perf.pair<-TableS7_Perf.pair %>% dplyr::rename(Timepoint = TimeP, Contrast =contrast, Estimate = estimate, DF = df, t = t.ratio, p = p.value, EffectSize = effect.size, EffectSize.SE = SE.es)
TableS7_Perf.pair<-TableS7_Perf.pair[,c("Timepoint", "Response", "Site", "Contrast", "Estimate", "SE", "DF", "t", "p", "EffectSize", "EffectSize.SE")]
#Round to 3 digits
TableS7_Perf.pair$Estimate<-round(TableS7_Perf.pair$Estimate, 3)
TableS7_Perf.pair$SE<-round(TableS7_Perf.pair$SE, 3)
TableS7_Perf.pair$t<-round(TableS7_Perf.pair$t, 3)
TableS7_Perf.pair$p<-round(TableS7_Perf.pair$p, 3)
TableS7_Perf.pair$EffectSize<-round(TableS7_Perf.pair$EffectSize, 3)
TableS7_Perf.pair$EffectSize.SE<-round(TableS7_Perf.pair$EffectSize.SE, 3)
#Integer
TableS7_Perf.pair$DF<-round(TableS7_Perf.pair$DF, 0)
##Write Out Table
write.csv(TableS7_Perf.pair, "Tables/TableS7_Performance_Pairwise_Results.csv", row.names=FALSE)
##Combine Results Tables
TableS8_IN_Tolerance<-data.frame(rbind(PAM_IN.res, Sym_IN.res, Chl_IN.res))
##Organize
names(TableS8_IN_Tolerance)
[1] "estimate" "SE" "df" "t.ratio" "p.value" "Predictor" "Response"
[8] "EtaSq"
TableS8_IN_Tolerance<-TableS8_IN_Tolerance %>% dplyr::rename(Estimate = estimate, DF = df, t = t.ratio, p = p.value, EffectSize = EtaSq)
TableS8_IN_Tolerance<-TableS8_IN_Tolerance[,c("Response", "Predictor", "Estimate", "SE", "DF", "t", "p", "EffectSize")]
#Round to 3 digits
TableS8_IN_Tolerance$Estimate<-round(TableS8_IN_Tolerance$Estimate, 3)
TableS8_IN_Tolerance$SE<-round(TableS8_IN_Tolerance$SE, 3)
TableS8_IN_Tolerance$t<-round(TableS8_IN_Tolerance$t, 3)
TableS8_IN_Tolerance$p<-round(TableS8_IN_Tolerance$p, 3)
TableS8_IN_Tolerance$EffectSize<-round(TableS8_IN_Tolerance$EffectSize, 3)
#Integer
TableS8_IN_Tolerance$DF<-round(TableS8_IN_Tolerance$DF, 0)
##Write Out Table
write.csv(TableS8_IN_Tolerance, "Tables/TableS8_Initial_Tolerance_LM_Results.csv", row.names=FALSE)